## **SpringerBriefs in Earth System Sciences**

**Olaf Kolditz · Keita Yoshioka · Tuanny Cajuhi · Ralf-Michael Günther · Holger Steeb · Frank Wuttke · Thomas Nagel Editors**

**GeomInt— Discontinuities in Geosystems From Lab to Field Scale**

## **SpringerBriefs in Earth System Sciences**

#### **Series Editors**

Gerrit Lohmann, Universität Bremen, Bremen, Germany

Justus Notholt, Institute of Environmental Physics, University of Bremen, Bremen, Bremen, Germany

Jorge Rabassa, Labaratorio de Geomorfología y Cuaternar, CADIC-CONICET, Ushuaia, Tierra del Fuego, Argentina

Vikram Unnithan, Department of Earth and Space Sciences, Jacobs University Bremen, Bremen, Bremen, Germany

SpringerBriefs in Earth System Sciences present concise summaries of cutting-edge research and practical applications. The series focuses on interdisciplinary research linking the lithosphere, atmosphere, biosphere, cryosphere, and hydrosphere building the system earth. It publishes peer-reviewed monographs under the editorial supervision of an international advisory board with the aim to publish 8 to 12 weeks after acceptance. Featuring compact volumes of 50 to 125 pages (approx. 20,000—70,000 words), the series covers a range of content from professional to academic such as:


Briefs are published as part of Springer's eBook collection, with millions of users worldwide. In addition, Briefs are available for individual print and electronic purchase. Briefs are characterized by fast, global electronic dissemination, standard publishing contracts, easy-to-use manuscript preparation and formatting guidelines, and expedited production schedules.

Both solicited and unsolicited manuscripts are considered for publication in this series.

Olaf Kolditz · Keita Yoshioka · Tuanny Cajuhi · Ralf-Michael Günther · Holger Steeb · Frank Wuttke · Thomas Nagel Editors

# GeomInt—Discontinuities in Geosystems From Lab to Field Scale

*Editors* Olaf Kolditz Environmental Informatics Helmholtz Centre for Environmental Research UFZ/TU Dresden Leipzig, Sachsen, Germany

Tuanny Cajuhi Underground Space for Storage and Economic Use Federal Institute for Geosciences and Natural Resources (BGR) Hannover, Lower Saxony, Germany

Holger Steeb Institute of Applied Mechanics University of Stuttgart Stuttgart, Baden-Württemberg, Germany

Thomas Nagel Geotechnical Institute TU Bergakademie Freiberg Freiberg, Sachsen, Germany Keita Yoshioka Environmental Informatics Helmholtz Centre for Environmental Research UFZ Leipzig, Sachsen, Germany

University of Leoben Leoben, Austria

Ralf-Michael Günther Institut für Gebirgsmechanik GmbH Leipzig, Sachsen, Germany

Frank Wuttke Geomechanics and Geotechnics Christian-Albrechts-University of Kiel Kiel, Schleswig-Holstein, Germany

ISSN 2191-589X ISSN 2191-5903 (electronic) SpringerBriefs in Earth System Sciences ISBN 978-3-031-26492-4 ISBN 978-3-031-26493-1 (eBook) https://doi.org/10.1007/978-3-031-26493-1

© The Editor(s) (if applicable) and The Author(s) 2023. This book is an open access publication. **Open Access** This book is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this book are included in the book's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the book's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use.

The publisher, the authors, and the editors are safe to assume that the advice and information in this book are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give a warranty, expressed or implied, with respect to the material contained herein or for any errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This Springer imprint is published by the registered company Springer Nature Switzerland AG The registered company address is: Gewerbestrasse 11, 6330 Cham, Switzerland

### **Acknowledgements**

The project "Geomechanical integrity of host and barrier rocks—experiment, modeling and analysis of discontinuities—GeomInt" is funded by the Federal Ministry of Education and Research (BMBF) within the geoscientific research funding programme "Geo:N Geosciences for Sustainability".

"Geo:N has an open structure. This facilitates identification of current issues over several years and redefinition of the respective priorities. The programme is intended to create stronger links between fundamental and applied geoscientific research and promote interdisciplinary research, Geo:N is part of the BMBF's Research for Sustainable Development (FONA)". More information can be found at www.ptj. de/en/project-funding/geo-n.

The funding under grant numbers 03G0866A-E is greatly acknowledged.


### **Contents**


### **Contributors**

**Markus Barsch** TUBAF, Technische Universität Bergakademie Freiberg, Freiberg, Germany

**Lars Bilke** UFZ, Helmholtz Centre for Environmental Research, OGS Core Team, Leipzig, Germany

**Jörg Buchwald** UFZ, Helmholtz Centre for Environmental Research, OGS Core Team, Leipzig, Germany

**Tuanny Cajuhi** BGR, Federal Institute for Geosciences and Natural Resources, Hanover, Germany

**Thomas Fischer** UFZ, Helmholtz Centre for Environmental Research, OGS Core Team, Leipzig, Germany

**Nico Graebling** UFZ, Helmholtz Centre for Environmental Research, Leipzig, Germany;

Leipzig University, Leipzig, Germany

**Norbert Grunwald** UFZ, Helmholtz Centre for Environmental Research, DE/OGS Core Team, Leipzig, Germany

**Carlos Guevara Morel** BGR, Federal Institute for Geosciences and Natural Resources, Hanover, Germany

**Ralf-Michael Günther** IfG, Institut für Gebirgsmechanik GmbH, Leipzig, Germany

**Nima Haghighat** CAU, Christian-Albrechts-Universität zu Kiel, Kiel, Germany

**Markus Knauth** IfG, Institut für Gebirgsmechanik GmbH, Leipzig, Germany

**Olaf Kolditz** UFZ, Helmholtz Centre for Environmental Research, Technische Universität Dresden, Dresden, Germany

**Dongwon Lee** UoS, University of Stuttgart, Stuttgart, Germany

**Christoph Lehmann** UFZ, Helmholtz Centre for Environmental Research, OGS Core Team, Leipzig, Germany

**Jobst Maßmann** BGR, Federal Institute for Geosciences and Natural Resources, Hanover, Germany

**Tobias Meisel** UFZ, Helmholtz Centre for Environmental Research, OGS Core Team, Leipzig, Germany

**Mostafa Mollaali** UFZ, Helmholtz Centre for Environmental Research, Leipzig, Germany

**Thomas Nagel** TUBAF, Technische Universität Bergakademie Freiberg, Freiberg, Germany

**Dmitri Naumov** UFZ, Helmholtz Centre for Environmental Research, OGS Core Team, Leipzig, Germany;

TUBAF, Technische Universität Bergakademie Freiberg, Freiberg, Germany

**Karsten Rink** UFZ, Helmholtz Centre for Environmental Research, OGS Core Team, Leipzig, Germany

**Amir S. Sattari** CAU, Christian-Albrechts-Universität zu Kiel, Kiel, Germany

**Patrick Schmidt** UoS, University of Stuttgart, Stuttgart, Germany

**Holger Steeb** UoS, University of Stuttgart, Stuttgart, Germany

**Wenqing Wang** UFZ, Helmholtz Centre for Environmental Research, OGS Core Team, Leipzig, Germany

**Frank Wuttke** CAU, Christian-Albrechts-Universität zu Kiel, Kiel, Germany

**Keita Yoshioka** UFZ, Helmholtz Centre for Environmental Research, Leipzig, Germany;

UoL, University of Leoben, Leoben, Austria

**Vahid Ziaei-Rad** UFZ, Helmholtz Centre for Environmental Research, Leipzig, Germany

**Gesa Ziefle** BGR, Federal Institute for Geosciences and Natural Resources, Hanover, Germany

### **Chapter 1 Introduction to GeomInt2**

#### **Thomas Nagel, Tuanny Cajuhi, Ralf-Michael Günther, Jobst Maßmann, Frank Wuttke, Keita Yoshioka, and Olaf Kolditz**

#### **1.1 Project**

GeomInt2 is the follow-up project to the joint project "Geomechanical Integrity of Host and Barrier Rocks—Experiment, Modelling and Analysis of Discontinuities (GeomInt)". The results of GeomInt made essential contributions to the analysis of the origin and development of discontinuities in clay, salt and crystalline rocks.

T. Nagel (B)

T. Cajuhi · J. Maßmann BGR, Federal Institute for Geosciences and Natural Resources, Hanover, Germany e-mail: Tuanny.Cajuhi@bgr.de

J. Maßmann e-mail: Jobst.Massmann@bgr.de

R.-M. Günther IfG, Institut für Gebirgsmechanik GmbH, Leipzig, Germany e-mail: ralf.guenther@ifg-leipzig.de

F. Wuttke CAU, Christian-Albrechts-Universität zu Kiel, Kiel, Germany e-mail: frank.wuttke@ifg.uni-kiel.de

K. Yoshioka UFZ, Helmholtz Centre for Environmental Research, Leipzig, Germany e-mail: keita.yoshioka@ufz.de

UoL, University of Leoben, Leoben, Austria

O. Kolditz UFZ, Helmholtz Centre for Environmental Research, Technische Universität Dresden, Dresden, Germany e-mail: olaf.kolditz@ufz.de

© The Author(s) 2023 O. Kolditz et al. (eds.), *GeomInt—Discontinuities in Geosystems From Lab to Field Scale*, SpringerBriefs in Earth System Sciences, https://doi.org/10.1007/978-3-031-26493-1\_1

TUBAF, Technische Universität Bergakademie Freiberg, Freiberg, Germany e-mail: thomas.nagel@ifgt.tu-freiberg.de

Discontinuities are considered in many forms, be it as volumetrically distributed damage macroscopically representing fissures in the rock microstructure, discontinuities that can form anew at phase interfaces as well as discrete cracks or fracture networks. The pathways created or extended by these discontinuities entail the risk of migration of fluid phases from deep to near-surface geological layers and must therefore be taken into account in geotechnical safety analyses.

The main result of the primarily methodology-oriented joint project GeomInt is the development of a broad experimental and numerical platform for laboratory and field analysis, and prediction of discontinuities in important reservoir and barrier rocks. Model comparisons for damage and fracture processes driven by different physical phenomena provide indications for the optimal application areas of the numerical methods. The follow-up project GeomInt2 aims to demonstrate the practical applicability of the developed experimental-numerical instruments under real conditions at the site scale. In addition, knowledge gaps that still exist will be closed through indepth tests on a smaller scale and in-depth model developments for specific process components.

At the core of GeomInt2 was an interdisciplinary consortium of partners from universities, state and private research institutions with complementary, long-standing experience in the analysis of geosystems that has proven itself in the previous project. New insights were elaborated, especially in understanding the effects of discontinuities on underground geosystems. In this context, the focus of experimental investigations shifted significantly to the implementation and evaluation of in-situ experiments in various underground research laboratories (URLs). For the efficient numerical simulation of coupled mechanical, thermal and hydraulic processes in the formation and development of discontinuities on the scale under consideration, the models and algorithms were further refined in GeomInt2. To enable the analysis of larger and more complex systems, high-performance computing (HPC) methods were added as a new aspect. The descriptive presentation of structural, experimental and model results in a real context (e.g. URLs) is to take place within the framework of an integrated visual data analysis approach (virtualisation). Figure 1.1 is a graphic illustration of the GeomInt2 project.

The project results aim at improved understanding of the processes, the used methods and the application-oriented systems for realistic time and length scales in order to make the planning and realisation of geotechnical uses of the subsurface safer, more reliable and more efficient. Another important advantage of GeomInt2 is the

#### **Fig. 1.2** GeomInt book

transferability of the experimental-numerical concepts and methods to other geotechnical applications such as deep geothermal energy, energy storage, repository problems, methods for hydraulic stimulation, conventional and unconventional resource extraction or tunnelling. GeomInt2 intensified the internationalisation started in the previous project by cooperation with complementary research projects (e.g. DECO-VALEX 2023 and the Mont Terri project) (Fig. 1.2).

The work packages mentioned in Fig. 1.1 were closely linked to the URLs Mont Terri (WP1/WP2), Springen (WP2) and Reiche Zeche (WP3) in continuation of the first project phase. The sample material for the geotechnical laboratories in Freiberg (TUBAF), Kiel (CAU) and Leipzig (IfG) originated from these URLs. Data from insitu experiments in Mont Terri (BGR) and the Springen pit (IfG) were incorporated into the modelling platform of the GeomInt project coordinated by the UFZ, in which all project partners participated. So-called Model EXperiments (MEX) were defined as a synthesis tool, in which experimental and modelling work were systematically combined from the beginning of phase one (Kolditz et al. 2021a, b). This integration formed the basis for the continuation of the project and the application of its results to site-scale problems in phase 2.

#### **1.2 Approach**

The research approach of GeomInt2 continues following the general concept of discontinuities in geosystems as illustrated in Fig. 1.3. Discontinuities in geological media—being a general feature in various host rocks such as salt, clay and crystalline rocks—occur at various scales. New fractures can be initiated by elevated fluid pressures, temperatures, and/or mechanical stresses in conjunction with geotechnical applications (TH2MCB drivers Fig. 1.3). Existing joints or faults can be reactivated

**Fig. 1.3** Graphical abstract of the GeomInt project: geological reservoir-barrier systems and geomechanical integrity (Nagel et al. 2021)

by similar multi-physical reasons as well (TH2MCB evolution Fig. 1.3). Fractures are localized and subject to size effects—they thus render a problem inherently multi-scale. Changes in connectivity of existing fracture networks due to evolving discontinuities can qualitatively alter system behaviour and lead to accelerated pathways from repositories to cap or side rocks (near-to-far-field bridging in Fig. 1.3).

*GeomInt2 leveraged both the experimental and numerical framework from the Lab to Field Scale*—which became the subtitle of the present book. GeomInt2 was aimed at improving our understanding of the underlying physical mechanisms and provide methods across different geological settings. The numerical instrumentation has been further developed and tested against experimental results on the field scale in particular from the underground research labs and mines in Mont Terri (CD-A experiment, Opalinus clay), Reiche Zeche (crystalline rock, STIMTEC experiment), and Springen (salt rock, large borehole experiment).

As a result of the GeomInt project so far, a broad experimental and numerical platform has been developed for the investigation of discontinuities due to swelling and shrinking processes (WP1), pressure-driven percolation (WP2) and stress redistribution (WP3). With these phenomena, processes relevant in the most important barrier and reservoir rocks (clay, salt, crystalline) are addressed. A comprehensive validation of the platforms ("Proof-of-Concept") was done by the Model EXperiments (MEX) for the damage and fracture processes driven by different processes. In the continuation of the project, the basic idea—the investigation of the three relevant process types for the formation of pathways in different rock types—is preserved. The proven project structure (WP1-3) is also retained (Fig. 1.4a).

**Fig. 1.4** GeomInt project geography

The emphasis in GeomInt2 is set as follows:


Figure 1.4 depicts the project geography of GeomInt2 in direct continuation of the workflows established in phase one of the project between the project partners from the universities of Freiberg (TUBAF), Kiel (CAU) and Stuttgart (UoS) as well as research institutions BGR, UFZ and the IfG company. The close collaboration established between the GeomInt and STIMTEC projects both working on discontinuities in crystalline rocks linked to the URL "Reiche Zeche" was also continued in the second phase.

The following chapters introduce some important results from the three work packages at the core of the project: WP1 (swelling/shrinkage), WP2 (fluid percolation) and WP3 (stress redistribution) for three types of host rocks, namely, clay, salt and crystalline rocks. Subsequently an introduction to the virtual reality and computation platforms is given, featuring the visualization, HPC and software engineering aspects. Some results are given more emphasis than others. Particularly, simulations on the hydro-mechanical behaviour of Opalinus clay were presented in greater detail. In the other chapters, a shorter presentation was chosen to keep the book concise. However, we point the reader to publications where more information also on the remaining topics addressed in the project can be found. The book closes with a synthesis and outlook.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 2 Hydro-Mechanical Effects and Cracking in Opalinus Clay**

**Tuanny Cajuhi , Nima Haghighat, Jobst Maßmann, Mostafa Mollaali, Amir S. Sattari, Vahid Ziaei-Rad, Gesa Ziefle, Thomas Nagel, Frank Wuttke, and Keita Yoshioka**

#### **2.1 Objectives and Outline**

In this chapter, we investigate hydro-mechanical effects in the Opalinus Clay, especially those leading to cracking. We present a methodology comprising laboratory and field scale experiments, as well as the development and application of numerical

T. Cajuhi (B) · J. Maßmann · G. Ziefle

J. Maßmann e-mail: jobst.massmann@bgr.de

G. Ziefle e-mail: gesa.ziefle@bgr.de

N. Haghighat · A. S. Sattari · F. Wuttke CAU, Christian-Albrechts-Universität zu Kiel, Kiel, Germany e-mail: nima.haghighat@ifg.uni-kiel.de

A. S. Sattari e-mail: amir.shoarian-sattari@ifg.uni-kiel.de

F. Wuttke e-mail: frank.wuttke@ifg.uni-kiel.de

M. Mollaali · V. Ziaei-Rad · K. Yoshioka UFZ, Helmholtz Centre for Environmental Research, Leipzig, Germany e-mail: mostafa.mollaali@ufz.de

V. Ziaei-Rad e-mail: vahid.ziaei-rad@ufz.de

K. Yoshioka e-mail: keita.yoshioka@ufz.de

K. Yoshioka UoL, University of Leoben, Leoben, Austria

© The Author(s) 2023 O. Kolditz et al. (eds.), *GeomInt—Discontinuities in Geosystems From Lab to Field Scale*, SpringerBriefs in Earth System Sciences, https://doi.org/10.1007/978-3-031-26493-1\_2

BGR, Federal Institute for Geosciences and Natural Resources, Hanover, Germany e-mail: tuanny.cajuhi@bgr.de

T. Nagel TUBAF, Technische Universität Bergakademie Freiberg, Freiberg, Germany e-mail: thomas.nagel@ifgt.tu-freiberg.de

approaches. In Sect. 2.2, we present a short overview concerning the Opalinus Clay (OPA) formation, its behavior and applications. This potential host rock for radioactive waste disposal shows complex coupled behavior, especially related to shrinkage and shrinkage/desiccation-induced cracking.

Section 2.3 describes the Cyclic-Deformation (CD-A) experiment in the Mont Terri Rock laboratory, one of the focal points of the GeomInt2 project. The experiment allows the measurement of relative humidity, deformation and cracking at in-situ scale. Furthermore, a drilling campaign allowed to extract core material that has been used in different laboratory experiments described in this chapter. We monitor the fracture response of OPA and evaluate the influence of its anisotropy (bedding orientation) as well as the effect of different suctions on the load-displacement curves. These experiments are modeled within numerical approaches based on the methods described in Sect. 2.4. In the latter, a continuum based, unsaturated hydro-mechanical approach is extended with the phase-field method (Sect. 2.4.2) and with the discrete element method (Sect. 2.4.3).

After comparing the response of the material both at the laboratory and at the field scales, we propose an interactive visualization strategy in Sect. 2.6. Finally, we review the proposed methodology in Sect. 2.7 and discuss its extension, future applications and transfer.

#### **2.2 Observations and Processes in Opalinus Clay**

The Opalinus Clay formation (OPA) is considered as a potential host rock for radioactive waste disposal due to its low hydraulic permeability and ability to retain radionuclides. Opalinus Clay formations are located in Northern Switzerland (Nagra 2021) and South-Eastern Germany (Hoth et al. 2007; Bundesgesellschaft für Endlagerung (BGE) 2020). In comparison to other clay formations, the OPA is well characterized by numerous in-situ experiments, conducted in the framework of the international Mont Terri Rock Laboratory in Switzerland. Its applicability in various geotechnical applications including the disposal of radioactive waste have been investigated, e.g. Bossart et al. (2018) and references therein.

The formation was accumulated during the dogger period (Jura) sandy and sandy facies as well as related subfacies types (Lauper et al. 2018; Kneuker et al. 2021). OPA is an anisotropic and highly heterogeneous material and shows complex and coupled thermo-hydromechanical (THM) behavior, including shrinkage and swelling. Swelling is related to an expansion of the material. In the current study, we focus on shrinkage, which is related to bulk volume reduction induced by environmental changes such as variation of relative humidity (RH), that, in turn is related to changes in water content. This effect can be observed and measured both in laboratory and field scales (Wild et al. 2017; Bock 2000; Zhang et al. 2010). Furthermore, it can lead to shrinkage-induced or desiccation cracks. With focus on how such cracks develop, the behavior of OPA is investigated in the following via experimental and numerical approaches that are applied at different scales. In-situ monitoring and borehole extraction are carried out in the Mont Terri Rock Laboratory in order to characterize the THM-behaviour of the OPA (Bossart et al. 2018; Kneuker and Furche 2021). In particular, the CD (Cyclic Deformation) experiment was a valuable starting point to gather experimental data and experience at the in-situ scale and to investigate the hydraulic-mechanical coupling induced by swelling and shrinking of shaly facies of the OPA due to cyclic variations of air humidity (Matray and Möri 2012; Ziefle et al. 2017). In the follow-up project, the CD-A-experiment, these processes were investigated further in detail, this time in the sandy facies of the OPA. Consequently, the CD-A-experiment is one focal point of the GeomInt2 project and will be described in the following.

#### **2.3 Monitoring and Data Integration**

In the Cyclic-Deformation (CD-A) experiment (Ziefle et al. 2022c), coupled effects due to the excavation and seasonal changes in relative humidity (RH) are monitored regularly. The long-term measuring program focuses on a detailed characterization of the host rock as well as on the evolution of process variables like the pore pressure, the displacements and the water content amongst others. The evolution of the excavation damaged zone (EDZ), corresponding changes of the water content and monitoring of desiccation cracks within the first decimeters of the rock are the central points of the experiment. These measurements enable an increased understanding of coupled HM effects in the EDZ and provide a valuable basis for a multi-disciplinary investigation aiming at further development of numerical modelling approaches and constitutive equations to reliability capture coupled effects in the OPA.

The CD-A experiment is characterized by two horizontal niches (twins) that are excavated perpendicular to the strike of bedding and that are not covered with any kind of stabilization (like shotcrete). Consequently, a direct interaction between the climatic conditions and the host rock takes place. While the HM effects in the open twin are evoked by the excavation and the seasonal variation of temperature and relative humidity, the relative humidity in the closed twin remains constant due to a sealing at its entrance. With that, the desaturation effect is minimized as much as possible holding the relative humidity high. The different evolution of the measured relative humidity within the niches is pictured in Fig. 2.1. The sudden drops in the red curve (closed niche) are related to door openings and installations. The open niche shows shows a cyclic decrease and increase of RH correspondent to the winter and summer months, respectively.

The changes in RH, especially those leading to an increase of suction, i.e., decrease in the capillary pressure, can lead to cracking. This phenomenon has been observed in both twins, but more intensely in the open niche. A detailed monitoring of the evolution of desiccation cracks along a horizontal line (scanline) in each niche is presented in Regard et al. (2021). The program includes monthly monitoring of the length and the aperture of the cracks initiated at the niche walls. The measurements have indicated significantly more desiccation cracks in the open than in the closed

**Fig. 2.1** CD-A experiment: measured relative humidity [%] in the twin niches. For example, the period between September 2019 and February 2020 shows a decrease of the RH, i.e. desaturation. Note that the RH in the closed niche remains high and approximately constant

twin. The sum of crack apertures in the open niche is also larger than in its closed twin. Due to the increase of relative humidity in the summer months, the sum of the crack apertures decreases.

In addition to the RH and cracking measurements, the CD-A experimental campaign comprises measurements focusing on the characterization of the claystone as well as the observation of process variables like the convergence, pore pressure and permeability which are amongst others measured in the boreholes presented in Fig. 2.2. Furthermore, laboratory scale tests are part of the study and provide further understanding on the OPA behavior. These tests have been carried out at material from CD-A drill cores extracted in March 2021 in the context of the permeability measuring campaign as presented in Fig. 2.2. Multi-disciplinary investigations along the presented boreholes are discussed in Ziefle et al. (2022a) and underline the high heterogeneity of the host rock with different facies and sub-facies types and the impact of the excavation process and the ventilation conditions on the near field around the niches.

In the following, we describe the performed laboratory experiments that focus on anisotropy characterization and mechanical response of the material under different suctions.

#### *2.3.1 Laboratory Experiments*

#### **2.3.1.1 Mechanical Characterization of Opalinus Clay**

In the first phase of the GeomInt project, fracture propagation characteristics in Opalinus clay have been investigated both numerically and experimentally. In this regard, a set of three-point-bending experiments has been performed on rectangular

**Fig. 2.2** Summary of boreholes from measuring campaign in 2022. Extensometers, minipiezometers and permeability sensors are installed, among others

samples of Opalinus clay. Given the status of bedding structure in obtained cores, only the configuration of loading being perpendicular to the bedding direction could be studied. Accordingly, the effect of structural anisotropy on sample strength as well as the fracture propagation behavior was investigated, albeit qualitatively, adopting different numerical methods. As an extension to the first phase of the project, a series of semi-circular bending experiments has been conducted in order to thoroughly investigate the direction-dependent strength characteristics of the Opalinus clay. The semi-circular configuration (Fig. 2.3) was chosen because of its minimal material requirements and relatively easier preparation procedure. The samples were taken from two different drill cores of sandy facies type, namely A24 and A26, near the CD-A twins (Kneuker and Furche 2021; Ziefle et al. 2022a, c).

Following ISRM suggested regulations (Kuruppu et al. 2014), the samples were prepared with radius of 36 × 10−<sup>3</sup> m, thickness of 40 × 10−<sup>3</sup> m and an artificial notch at center with length of 14 × 10−<sup>3</sup> m. Given sensitivity of the clay material to water, all specimens were prepared under dry conditions using a diamond band saw with 2 × 10−<sup>3</sup> m thickness. As summarized in Fig. 2.3, three type of samples, namely P-, Z- and S-samples, have been considered for the experiments. The naming convention is indicative of loading angle θ, which represents the orientation between the isotropy plane and the loading direction. We have tested following loading directions θ = [0, 45, 90] ◦ that represent P-, Z- and S-samples, respectively.

In total, three samples, two from A26 cores and one from A24 core, for each loading angle were tested. The tests were performed under quasi-static loading conditions with strain rate of 0.03 mm/min. The typical laboratory results in terms of applied load versus the deflection of the samples obtained from A24 core drill are illustrated in Fig. 2.4a. The load-deflection curves in the three cases exhibited a similar pattern; a downward convex trend followed by an approximately linear increase

**Fig. 2.3** Schematic summary of SCB samples: (from left to right) P-samples, Z-samples and S-samples

**Fig. 2.4 a** Obtained Load-Deflection curve for Opalinus clay sample with different loading angles, **b** Mean value and range of the obtained peak loads for each loading angle. With increasing loading angle, an increase of the peak load is observed

in the load response until failure. Furthermore, the sharp drop of the load represents the brittle failure mechanism during the experiments. The mean value and the range between the maximum and minimum achieved failure loads are summarized in Fig. 2.4b. As shown, the maximum obtained load before failure, indicative of material toughness, increases as the angle between the loading direction and the isotropy plane increases. The considerable range among toughness of the samples with similar anisotropy structure could be attributed to various factors, e.g. micro defects generated in samples during the preparation process or the conditions at which cores were drilled.

Figure 2.5 shows the final fracture pattern of the three samples from A24 cores. Note that the crack does not propagate straight along the bedding direction for the P-sample due to the heterogeneity of the material. In the Z-sample, the crack propagation along bedding direction can be identified more clearly. Finally, the crack propagation in the S-sample also shows, initially, a strong dependency on the bedding direction, but continues its trajectory straight, as expected.

**Fig. 2.5** Failure pattern and fracture path trajectory of **a** P-sample, **b** Z-sample and **c** S-sample. The cracks propagate approximately straight from the notch

#### **2.3.1.2 OPA Characterization Under Different Humidity Conditions**

As also observed at in-situ scale, in the CD-A experiment, changes in the RH lead to variations in the mechanical and fracture behaviors. We investigate this phenomenon in the laboratory by characterizing the mechanical response of the OPA at various RH. To do that, a customized desiccator cell has been configured and set up at the geomechanics and geotechnics laboratory of CAU. The idea originates from the first phase of the project, where the concept of a humidity controlled bending test setup was proposed (Sattari et al. 2021). As shown in Fig. 2.6, the customized desiccator cell contains four distinct parts, namely


The test procedure generally takes place in two steps; the equilibrium phase and the loading step. The sample is exposed to a certain RH during the equilibrium phase while being allowed to adjust with the environment. The humidity conditions inside the cell were achieved and controlled by applying different supersaturated salt solutions. For the Opalinus clay samples, the equilibrium normally occurs between 7 to 10 days. The equilibrium period varies due to different initial water content of the sample and the specified RH of the cell. The equilibrium conditions are monitored with consecutive weight measurements of each sample. The unique design of the cell provides an opportunity to assess mechanical characteristics of OPA under controlled humidity conditions with minimizing the environmental exposure of sample between equilibrium phase and loading step. In this section, anisotropic strength characteristics of Opalinus clay are studied in different humidity conditions using the semi-circular bending configuration. Table 2.1 summarizes the four different humidity conditions considered in this study. The obtained suction pressure for each

**Fig. 2.6** Customized Desiccator Cell equipped with relative humidity (RH) and temperature sensors, a salt bath, an interchangeable mechanical loading apparatus and an adjustable camera tube


**Table 2.1** Utilized salt solutions along with the achieved relative humidity and generated suction

salt solution using 3 samples each is estimated using Kelvin's equation and based on the recorded RH and temperature during the experiments,

$$\psi = -\frac{R \, T}{V\_{w0} \, M} \ln \left( \frac{p}{p\_0} \right) \tag{2.1}$$

where the universal gas constant *R* = 8.314 J/(mol K), *T* is the temperature in Kelvin, *V*w<sup>0</sup> is the specific volume of water and *M* = 0.01801528 kg/mol is the molar mass of the pore fluid. The term *p*/*p*<sup>0</sup> corresponds to the relative humidity (RH).

The obtained experimental results are presented in terms of failure load as a function of imposed suction for each loading angle in Fig. 2.7. A general increasing trend for the material strength following the increase in the imposed suction and resultant decrease of water content can be clearly observed. This behavior is associated with

**Fig. 2.7** Load-displacement curves under different RH. In general, the load increases with suction

the generation of capillary forces in a partially saturated porous medium due to the pressure difference between gas and fluid phases at the pores (Wild et al. 2015).

Considering the aspects mentioned in this section, we aim at reproducing the shrinkage-behavior of OPA numerically in the following.

#### **2.3.1.3 Additional Tests on Percolation**

We utilized the drilled cores to carry out fluid percolation tests. These tests are important in the context of hydraulic fracturing. Knowing when the fluid will induce a crack aids in understanding critical conditions in the field. Fluid percolation experiments were carried out in the geomechanics and geotechnics laboratory of CAU Kiel. The main objective to these experiments is to investigate the characteristics of fluid driven fractures in Opalinus clay under effects of structural anisotropy as well as the anisotropy of the stress state. The utilized setup includes a true triaxial loading apparatus, resembling the in-situ stress conditions, and a syringe pump. More details related to experimental setup and the true-triaxial device can be found in Sattari et al. (2021). Cubic samples were prepared with length of 43 × 10−<sup>3</sup> m and a scaled borehole drilled with diameter of 7 × 10−<sup>3</sup> m and length of 20 × 10−<sup>3</sup> m. Fluid is injected under constant flow rate into the sample through an adapter and a steel tube glued to the drilled borehole. A fluorescent oil was used as the injection fluid in order to trace even the tiniest fractures initiated during the process. Schematic representation

**Fig. 2.8** Schematic representation of cubic S-, P- and Z-samples (left to right)


**Table 2.2** Summary of fluid percolation experiments on Opalinus clay

of the applied principle stresses along with the bedding direction of the samples are shown in Fig. 2.8.

Given the position of the scaled borehole with respect to the bedding direction, three different loading configurations can be tested:


Accordingly, Table 2.2 summarizes the number of tested samples for each bedding direction and the stress state in which the experiment was conducted.

The experimental test starts after the in-situ stresses are applied. Once the mechanical equilibrium has been reached, the fluid injection then starts at flow rate of 10 ml/h. Upon observing the breakdown pressure, the experiment is stopped and the sample is retrieved. After careful investigation of failure patterns on the surface, the samples were sawed into number of cross section to have a better understanding of fracture propagation behavior inside it. Figure 2.9 shows the crack morphology on the surface, the outline of the cross sections, and the inner trajectories of fracture path indicated by florescent oil under UV light. It can be clearly seen that the pressure induced fractures inside Opalinus clay can either follow a straight path or propagate in branches and form a network. These observations are also reflected in the pressure response shown in Fig. 2.10. Generally, the drop in the pressure response is associated with increase in the occupied volume by the injection fluid followed by initiation and propagation of fractures inside the sample. In some cases, e.g., Fig. 2.10c, the initiated fracture reaches the surface of the sample rapidly and causes the depletion of the injection fluid out of the sample and a sharp breakdown in pressure response. In other cases, e.g., Fig. 2.10b, fractures propagate in tiny branches and form a network before reaching the sample surface which is reflected in the fluctuations of the pressure response.

2 Hydro-Mechanical Effects and Cracking in Opalinus Clay 17

**Fig. 2.9** Fracture path trajectory, cross-section cut and inner fracture path trajectory indicated by florescent oil for **a** P-sample **b** S-sample and **c** Z-sample of Opalinus clay under isotropic in-situ stress of 12MPa

**Fig. 2.10** Comparison between pressure-time curve for different isotropy planes of Opalinus clay under mechanical stress configurations of **a** *S*<sup>v</sup> = 12 MPa, *S*<sup>H</sup> = 12 MPa, *S*<sup>h</sup> = 12 MPa, **b** *S*<sup>v</sup> = 12 MPa, *S*<sup>H</sup> = 8 MPa, *S*<sup>h</sup> = 8 MPa, and **c** *S*<sup>v</sup> = 12 MPa, *S*<sup>H</sup> = 8 MPa, *S*<sup>h</sup> = 4 MPa

#### **2.4 Numerical Modeling Approach**

To reproduce shrinkage effects observed in OPA in the laboratory and field scales (CD-A experiment), we utilize and extend the numerical methods evaluated in the first phase of the GeomInt project (Yoshioka et al. 2021; Kolditz et al. 2021a). The extensions are related to the anisotropic and cracking behaviors of the material.

Intuitively, cracks generate edge separations with creation of boundary surfaces over which the displacements are discontinuous. This makes numerical simulation of crack propagation challenging especially if we do not prescribe the location, topology, or time of their creation. For modeling crack propagation, several strategies are proposed in the literature. One way to deal with this sort of topology changes is by adaptively remeshing the domain (Bouchard et al. 2000; Meyer et al. 2004; Branco et al. 2015; Rangarajan et al. 2015). However, assigning boundary conditions to the newly created surfaces and including potential contact are not straightforward tasks. Another way to model crack propagation is to enrich elements that contain cracks to capture the displacement discontinuity, an approach known as eXtended/Generalized Finite Element Methods (XFEM, GFEM) (Moës et al. 1999; Belytschko et al. 2001; Gasser and Holzapfel 2005; Fries and Belytschko 2006; Belytschko et al. 2009; Khoei et al. 2012). Despite recent success (Shauer and Duarte 2022), the implementation efforts remain a daunting task especially when arbitrary fracture topologies are considered (Belytschko et al. 2001, 2003; Duarte et al. 2007).

Crack initiation needs to be prescribed in the methods discussed above. This can become a limitation in applications of geomechanical integrity, where we need to assess fracture initiation possibilities without any a-priori knowledge. The phasefield approach for brittle fracture represents the crack using a phase-field (damage) variable that results in a diffuse cracked interface. Furthermore, it does not require any special mesh treatment and prescription of the crack position. We apply and evaluate the response of both diffuse cracks using a continuum mechanics based approach combined with the phase-field and sharp cracks using the discrete element method. While the former is based on successive minimization of the total energy and cracks initiate as a result of the energy minimization, the latter finds fracture initiation based on the stress as the stress exceeds the strength defined between discrete elements.

#### *2.4.1 Unsaturated Hydro-Mechanical Approach*

The continuum formulation of the unsaturated HM-coupled problem outlined below is discussed in detail in Pitz et al. (2022) and is mainly based on a simplification of the approach from Grunwald et al. (2022). It is implemented in the current version of the open source code OpenGeoSys (Bilke et al. 2022).

The isothermal formulation is founded on the conservation of water volume in a deformable porous medium according to Biot's consolidation theory (Biot 1941) considering unsaturated Darcy flow by the Richards equation (Richards 1931):

$$\begin{split} 0 &= \rho\_{\text{LR}} \text{S}\_{\text{L}} \left( \phi \beta\_{p,\text{LR}} + \text{S}\_{\text{L}} (\alpha\_{\text{B}} - \phi) \beta\_{p,\text{SR}} \right) \frac{\text{d} p\_{\text{LR}}}{\text{d} t} - \text{div} \left( \rho\_{\text{LR}} \frac{k\_{\text{L}}^{\text{eff}} \mathbf{k}\_{\text{S}}}{\mu\_{\text{LR}}} \left( \text{grad } p\_{\text{LR}} - \rho\_{\text{LR}} \mathbf{g} \right) \right) \\ &+ \left[ \rho\_{\text{LR}} \phi + \rho\_{\text{LR}} p\_{\text{LR}} \mathbf{S}\_{\text{L}} (\alpha\_{\text{B}} - \phi) \beta\_{p,\text{SR}} \right] \frac{\text{d} \mathbf{S}\_{\text{L}}}{\text{d} t} + \text{S}\_{\text{L}} \rho\_{\text{LR}} \alpha\_{\text{B}} \text{ div } \frac{\text{d} \mathbf{u}}{\text{d} t} \end{split} \tag{2.2}$$

with the density of the liquid phase ρLR, the liquid saturation *S*L, the porosity φ, the liquid and solid compressibility β*<sup>p</sup>*,LR, β*p*,SR, respectively, the Biot coefficient αB, the liquid pressure *p*LR, the intrinsic permeability tensor **k**<sup>S</sup> and relative permeability *k*rel <sup>L</sup> , the liquid viscosity μLR and the solid displacement vector **u**.

The capillary pressure *p*cap is defined as the difference between the liquid pressure *p*LR and the gas pressure *p*GR. Following the Richards equation, pressure gradients in the gas phase are neglected and by choosing the gas pressure to zero, we obtain following relationship

$$p\_{\rm cap} = p\_{\rm GR} - p\_{\rm LR} = -p\_{\rm LR} \,. \tag{2.3}$$

The liquid saturation *S*<sup>L</sup> is linked to the effective liquid saturation via *S*<sup>L</sup> = *S*e *S*L,max − *S*L,res + *S*L,res, which is in turn defined as a function of the capillary pressure according to van Genuchten (Van Genuchten 1980):

$$S\_{\mathfrak{e}} = \left( 1 + \left( \frac{p\_{\rm cap}}{p\_{\mathfrak{e}}} \right)^n \right)^{\frac{1}{n} - 1},\tag{2.4}$$

with *p*<sup>e</sup> as the gas entry pressure and *n* ∈ [0, 1) as empirical material property. The relative permeability is defined according to van Genuchten/Mualem (Mualem1976):

$$k\_{\rm L}^{\rm rel} = \sqrt{S\_{\rm e}} \left( 1 - \left( 1 - S\_{\rm e}^{n/(n-1)} \right) \right)^{(2n-2)/n} . \tag{2.5}$$

The momentum balance equation is given by

$$\mathbf{0} = \text{div}\left(\sigma' - \alpha\_{\text{B}} \mathbf{S}\_{\text{L}} p\_{\text{LR}} \mathbf{I}\right) + \left[\left(1 - \phi\right)\rho\_{\text{SR}} + \mathbf{S}\_{\text{L}} \phi \rho\_{\text{LR}}\right] \mathbf{g} \,, \tag{2.6}$$

where ρSR is the density of the solid phase, **g** the earth's gravitational acceleration vector and σ represents the effective stress tensor, which can be calculated by the solid stiffness **C** and the elastic strains: σ = **C**: ( − inelastic). In the current approach, inelastic strains are not taken into account, i.e. inelastic = 0.

#### *2.4.2 Extension with the Phase-Field Approach (Diffuse Cracks)*

The coupled HM formulation is extended to account for cracking with the phasefield modeling approach for brittle fracture. With the phase-field approach (Bourdin

**Fig. 2.11** Discrete representation of fractures (left) and diffused phase-field representation of fractures (right)

et al. 2000, 2008), we introduce a scalar variable that represents the damage of the material *d* : → [0, 1]. A material, which is free from damage (intact) is represented by *d* = 0, while a fully degraded material by *d* = 1. Damage evolves in a diffuse manner, i.e. the crack phase-field variable evolves continuously, as shown on the right side of Fig. 2.11.

We couple the evolution of the crack phase-field *d* with the momentum balance equation by multiplying a degradation function *g*(*d*) based on the evolution of *d* with the effective stresses. Before doing so, we split the effective stresses into active σ <sup>+</sup> and non-active σ <sup>−</sup> parts computed from decomposing the elastic strain energy (Amor et al. 2009; Miehe et al. 2010) as follows

$$\begin{split} \sigma' &= g(d) \frac{\partial \Psi\_+}{\partial \epsilon} + \frac{\partial \Psi\_-}{\partial \epsilon} \\ &= g(d) \sigma'\_+ + \sigma'\_- \end{split} \tag{2.7}$$

Note that only the active, crack-driving part is degraded by *g*(*d*). The momentum balance equation is rewritten as

$$\mathbf{0} = \operatorname{div} \left( \mathbf{g}(d) \sigma\_+^{\prime} + \sigma\_-^{\prime} - \alpha\_\text{B} \mathbf{S}\_\text{L} p\_\text{LR} \mathbf{I} \right) + \left[ (1 - \phi) \, \rho\_\text{SR} + \mathbf{S}\_\text{L} \phi \rho\_\text{LR} \right] \mathbf{g} \,. \tag{2.8}$$

In this study, *g*(*d*) is given by

$$g(d) = (1 - d)^2 + \eta,\tag{2.9}$$

where η is the residual stiffness (e.g., η = 1 × 10−<sup>8</sup> Pa) assigned in fully degraded part (*g*(*d*) = 0) for numerical stability. An overview regarding the derivation and application of a similarly combined unsaturated HM and phase-field approaches is described in Cajuhi et al. (2018). The crack phase-field evolution is given by

#### 2 Hydro-Mechanical Effects and Cracking in Opalinus Clay 21

$$g'(d)\Psi\_{+}(\mathbf{u}) + \frac{\mathcal{G}\_c}{4c\_w} \left(\frac{w'(d)}{\ell} + 2\ell\Delta d\right) = 0,\tag{2.10}$$

where

$$w(d) := (1 - d)^n \quad \text{and} \quad c\_w := \int\_0^1 (1 - s)^{n/2} \, \text{ds},\tag{2.11}$$

and *n* = 1 and 2 stand for so-called AT<sup>1</sup> and AT<sup>2</sup> models respectively (Pham et al. 2011; Bourdin et al. 2014), and *c*<sup>w</sup> = 2/3 for AT<sup>1</sup> and 1/2 for AT2.

In the phase-field literature, the split of the energy is sometimes referred to anisotropic response (Ambati et al. 2015). In this book, the term anisotropy refers to the bedding orientation and constitutive behavior. For referencing the energy decomposition (split), we use the term tension-compression asymmetry. There are a few popular split models in the literature, nevertheless, most of them can only be used for isotropic materials. As for anisotropic materials, the principal directions of the stress and strain no longer coincide. Hence, the orthogonality condition, a requirement for the split method to be variationally consistent, is not guaranteed. Thus, further treatment is required. In our recent work, we generalized two existing split models for isotropic materials in order to account for materials with anisotropic constitutive behavior. Borrowing some basic concepts from He and Shao (2019), we performed an energy preserving transformation so that the orthogonality condition is satisfied, and thus, a decomposition of an arbitrary anisotropic elastic model is feasible. We encourage those interested readers to follow our manuscript for a detailed clarification (Ziaei-Rad et al. 2023). Exemplarily, we show the verification of the framework in Sect. 2.5.1.2.

#### *2.4.3 Extension with the Discrete Element Method (Discrete Cracks)*

The coupled HM formulation is extended to account for cracking with the Finite Discrete Element Method (FDEM).

#### **2.4.3.1 FDEM Mechanical Solver**

Originally proposed by Munjiza (2004), the Finite Discrete Element Method (FDEM) utilizes principles of continuous numerical approaches in combination with those of discontinuous methods to simulate interaction of multiple deformable bodies (Lisjak et al. 2014). In other words, while the elastic deformation of discrete bodies and the initiation and propagation of fractures are captured using Finite Element analysis, the interaction between individual blocks of elements are realized by Discrete Element analysis. In FDEM, each block of the problem domain is discretized into a finite

**Fig. 2.12 a** FDEM discretization of the intact discrete blocks containing triangular finite elements and zero-thickness joint elements, **b** dual grid of FDEM flow solver

element mesh with zero-thickness joint elements being inserted on its intersection (common) edges of the adjacent elements as shown in Fig. 2.12. Then, an explicit finite difference time integration scheme is applied to solve the equation of motion for the discretized system which reads as

$$\mathbf{M}\frac{\partial^2 \mathbf{x}}{\partial t^2} + \mathbf{C}\frac{\partial \mathbf{x}}{\partial t} + \mathbf{F}\_{\text{int}}(\mathbf{x}) - \mathbf{F}\_{\text{ext}}(\mathbf{x}) - \mathbf{F}\_c(\mathbf{x}) = 0,\tag{2.12}$$

where**M** and **C** denote the diagonal lumped mass and damping matrices, respectively. **x** is the nodal displacement vector, and *F*int, *F*ext and *Fc* are the nodal internal, external and contact forces, respectively. Contact forces, *Fc*, are calculated between contacting discrete blocks either in the normal direction, using a penalty method, or in the tangential direction, employing Coulomb-type friction law. Accordingly, internal resisting forces, *F*int, are calculated as summation of elastic forces and joint elements' bonding forces.

Following a cohesive-zone approach, the mixed-mode crack model consists of a typical stress-strain curve where the bonding stresses are proportional to the separation and the sliding of the fracture edges. More details of the crack model can be found in the literature (Mahabadi et al. 2012).

#### **2.4.3.2 FDEM Flow Solver**

In order to account for the saturated/unsaturated porous media flow within the framework of FDEM, a vertex-centered finite volume scheme is adopted. Based on this approach, the conservation of the fluid phase is enforced on a dual grid of polygonal control volumes. As illustrated in Fig. 2.12, each control volume is generated around vertices of the mechanical solver grid connecting centers of triangular elements through midpoints of primary grid. The adopted mass balance equation of fluid phase in the porous medium is a simplified version of Eq. 2.2 with incompressible solid phase. The hydro-mechanical coupling in FDEM is realized by a staggered approach alternating between solid and flow solvers. Within this approach, the valuation of the stress field is affected by the pore pressure according to the effective stress principle, while the flow calculations are affected by the volume change associated with each computational node.

#### **2.5 Process Simulation Results**

We show different simulation setups and their results using the methods described in Sect. 2.4. At first, we model the mechanical processes at laboratory scale considering the anisotropy of the Opalinus Clay and using the FDEM formulation (Sect. 2.4.3). The response of the samples is compared with the laboratory experiments discussed in Sect. 2.3.1. Cracks are modeled either discretely or diffusely and we compare both approaches. At the field scale, the deformation and pore pressure behaviors are described with the continuum approach from Sect. 2.4.1 and the hydro-mechanical effects such as the shrinkage due to the de-saturation are discussed. We model the formation of the excavation damaged zone (EDZ) with discrete cracks (formulation from Sect. 2.4.3). Finally, desiccation cracks in the field are modeled with the combined poromechanical phase-field approach (formulation from Sect. 2.4.2).

#### *2.5.1 Laboratory Scale*

#### **2.5.1.1 Crack Paths Due to Deformation**

As indicated by the achieved experimental results in Sect. 2.3.1.1, not only the deformation response but also the strength characteristics of Opalinus clay are strongly influenced by the bedded structure of the samples. In the 2D FDEM method, the anisotropic elastic behavior is captured by a transversely isotropic stress-strain law, while the strength directionality is modeled adopting a so-called smeared approach (Lisjak et al. 2014). Within this approach, the strength of each joint element is directly linked to its relative orientation with respect to the bedding layers. In this way, the strength parameters and fracture energies fluctuate linearly between the minimum limiting value, when the cohesive element is parallel to the bedding layer, and the maximum limiting value, when the cohesive element is perpendicular to the bedding layer.

The FDEM models were calibrated to match the obtained experimental peak loads reported in Sect. 2.3.1.1. The laboratory-scale model along with the employed boundary conditions as well as the mesh details for P-, S- and Z-samples are illustrated in Fig. 2.13. The experimental loading conditions were realized by explicitly

**Fig. 2.13 a** Numerical model setup of SCB experiment and mesh details for **b** P-sample **c** Z-sample and **d** S-sample


**Table 2.3** Material and input parameters for FDEM simulations of three point bending test of OPA

modeling the supports and the loading platen. In order to avoid the impact effect of loading platen, the loading velocity increases linearly from zero to 0.1 m/s in 1 ms. Following the procedure reported in the work of Tatone and Grasselli (2015), the calibrated mechanical material properties in the FDEM model are summarized in Table 2.3.

**Fig. 2.14** Fracture path trajectory (top) and corresponding distribution of horizontal stress [Pa] (bottom) at 0.5 mm of crack mouth opening displacement for Opalinus clay of P-, Z- and S-samples with **a** θ = 0◦, **b** θ = 45◦ and **c** θ = 90◦, respectively

Figure 2.14 illustrates the numerical simulation results concerning the fracture path trajectories upon failure and the distribution of horizontal stress for the three cases. The crack propagation shows good agreement with the experimental observations depicted in Fig. 2.5.

#### **2.5.1.2 Pure Shear Test**

In order to verify the anisotropic phase-field model detailed in Ziaei-Rad et al. (2023), we present a set of numerical examples of a single edge notched shear test. The geometry and boundary conditions for this examples are shown in Fig. 2.15. We consider a square plate (0.001 m × 0.001 m) under plane strain conditions with a pre-notched crack at the middle height of the specimen. The length scale parameter was chosen as = 1 × 10−<sup>5</sup> m which is small enough with respect to the specimen dimensions. The AT<sup>1</sup> model was used in the simulations. We used the parameters of Opalinus clay listed in Table 2.3. Note that the critical surface energy *Gc* = 100 N/m was taken to be independent of the material orientation.

As boundary condition, we apply fixed displacement (zero in all directions) at the bottom edge (*y* = −*L*/2). For simulating shear, we prescribe a displacement along the *x*-direction at the top edge (*y* = *L*/2). The computation was performed with a monotonic displacement controlled loading with a constant displacement increment *u* = 3 × 10−<sup>7</sup> m.

Figure 2.16 depicts the crack paths for orthotropic volumetric-deviatoric (Amor et al. 2009) split under shear loading with three material orientation angles

**Fig. 2.15** Schematic of a cracked square plate under a displacement shear loading. An incremental displacement loading is applied on the top edge. The material principal axes (*e*1, *e*2) have an angle of α with the global coordinate system (*x*, *y*)

(α = −π/4, 0, and π/4). Herein a significant difference in crack response is observed between the three cases.

The volumetric-deviatoric decomposition assumes the volumetric compression does not contribute to the crack evolution. However, in case of anisotropy, the Poisson effect may pronounce more in a direction compared to the others, thus, even under a uniaxial compression, some crack driving forces can lead into damage development. Consequently, the crack response of a compression-dominated loading is depending on differences in directional stiffness of material. An implementation of the presented anisotropic phase-field formulation combined with the unsaturated HM approach is underway.

#### *2.5.2 Field Scale*

#### **2.5.2.1 Hydro-Mechanical Effects and Model Response**

In this section, we numerically model the hydro-mechanical coupled effects in the context of the CD-A experiment using the macroscopic approach presented in

**Fig. 2.16** Cracked square plate under shear. Phase field contours of orthotropic volumetricdeviatoric for three material orientations. The same parameters were used for all cases. Here, we observe that the crack paths are strongly dependent on the material orientation. With α = π/4 there exist straight cracks that grow towards the lower right corner, almost along the material principal axis with lower stiffness. The two other cases differ from an expected isotropic response

Sect. 2.4.1. The two-dimensional model consists of a cross section of the open and closed niches and surrounding domain, Fig. 2.17. We prepare an isotropic and an anisotropic model, in which the bedding direction is assumed to be approximately horizontal as monitored in the field (Ziefle et al. 2021). The OPA parameters and properties are extracted from several references in the literature (Bock 2000; Martin and Lanyon 2003; Garitte et al. 2017; Bossart et al. 2018; Zhang et al. 2019) and presented in Table 2.4. The transversal component of an anisotropic property is defined with the subscript 1, while its parallel components are defined with subscripts 2 and 3. These subscripts corresponds to a local coordinate system oriented according to the bedding layers, i.e. 1 is orthogonal to the bedding. We compute the isotropic equivalent properties as the average of the anisotropic properties, e.g. the isotropic equivalent elastic modulus corresponds to the average of its parallel and perpendicular components. The isotropic equivalent properties are denoted with the subscript (.)iso in Table 2.4. The initial stress values are based on in-situ stress measurements (Martin and Lanyon 2003). Their orientation correspond to the global coordinate system, i.e. in two-dimensional setting their horizontal and vertical components σ<sup>h</sup> and σ<sup>v</sup> lie on the global coordinate system (*x*, *y*) (Fig. 2.17), while in three-dimensional setting σh, σ<sup>H</sup> and σ<sup>v</sup> on (*x*, *y*,*z*) (Fig. 2.21). The initial pore pressure and in-situ stresses are applied on the domain. After instantaneous excavation, the seasonal boundary condition computed from the RH measurements (Fig. 2.1) with the Kelvin's equation (Eq. 2.1) is applied at the niches. In the closed niche, the drops in the RH curve are related to installations and openings of the sealing door. We consider that the latter remains under high RH, i.e. neglect the drops, and assume a constant pore pressure of −1 MPa. The open niche follows a cyclic pressure evolution, interpolated from the first season of the RH curve, Fig. 2.1.

Figure 2.18 shows the saturation results using the anisotropic model for March 2020 (first desaturation peak in winter) Fig. 2.18a and September 2020 (first re-


**Table 2.4** Parameters and material properties used for simulating the CD-A experiment

saturation peak in summer) Fig. 2.18b. Note that the horizontal direction is more permeable and, consequently, the de-saturated zone is more extended in this direction when compared to its orthogonal counterpart. While the saturation in the closed niche remains high, the saturation near the wall of the open niche increases between March and September. At the same time, the overall desaturated zone extends along the domain, but still indicates a local effect. Except for direction dependency, the same conclusions are drawn in the isotropic model.

The aforementioned tendencies, i.e. de-saturation and re-saturation during winter and summer, respectively, have been compared with the monitored water content

#### 2 Hydro-Mechanical Effects and Cracking in Opalinus Clay 29

**Fig. 2.17** Setup of the two-dimensional model of the CD-A experiment consisting of open and closed niches. The niche radius is 1.15 m and distance between them is 10 m. The bedding direction is horizontal and indicated in the figure

**Fig. 2.18** Saturation changes in deformation results in the 2D model of the CD-A experiment. The region near the open niche dries, de-saturates, during winter and re-saturates during the summer months. Convergence is observed due to an increase of the deformation near the niches. Note that both effects are more pronounced in the open niche in comparison to the closed niche

and electrical resistivity results. A decrease in saturation has been associated to a decrease of water content and an increase of the resistivity, while the opposite occurs when the saturation increases (Amabile et al. 2020; de Jong et al. 2020; Cajuhi et al. 2021). These comparisons, including a geologic characterization of the niches have revealed that inhomogeneities like cracks, shear zones, carbonate-rich and sandy lenses subfacies in the claystone have a strong influence on the seasonal desaturation or re-saturation behavior in the near field. A detailed description of this study and its conclusions are reported in Ziefle et al. (2022a, b).

The numerical results show that the seasonal saturation changes affect the deformation behavior. We compare the strain response of the isotropic and anisotropic models, Figs. 2.19 and 2.20, respectively, in the first winter and summer peaks. On both niches, the contribution of the strains due to excavation is equal. The increase of deformation represents convergence of the niches. Note that the projected initial stress is the same in both isotropic and anisotropic models. In the isotropic model, the strain responses are similar, while the strains in vertical direction are larger than in the horizontal direction in the anisotropic model. Note that the strains decrease from winter to summer when comparing the responses of both models. This decrease is specially observed in the vertical and shear components, and particularly, in the

**Fig. 2.19** Deformation response of isotropic model. Blue regions denote negative strains, while red ones relate to positive strains. The strain during the first desaturation peak (March 2020) are shown on the left, while the strains during the second peak are shown on the right (September 2020). An equal strain response is observed near the niches due to excavation. The change of saturation induces slightly larger deformation in the vicinity of the open niche

open niche, where seasonal effects are more pronounced and accounted for in the model.

As the driving mechanism of the crack phase-field is the active elastic energy, which evolves with the strains (see Sect. 2.4.1), we expect to model cracks due to drying in the open niche using this approach. Desiccation cracks are expected to occur due to mode I loading, i.e. tension cracks. In the field, they were observed on the wall surfaces, i.e. where the rock is exposed and consequently drier in comparison to the undisturbed rock. Such cracks are expected to propagate from the niche walls into the depth of the rock. To further understand the stress distribution on the wall surfaces, we have prepared a 3D model of the experiment. For completeness, the model also considers the lock to the closed niche and the gallery that provides access to the niches. The seasonal boundary conditions are applied in the open niche, lock and gallery. The numerical model is shown in Fig. 2.21. The numerical model has been also utilized for interactive visualization purposes in the virtual laboratory

**Fig. 2.20** Deformation response of anisotropic model. Blue regions denote negative strains, while red ones relate to positive strains. The strains due to the first desaturation peak (March 2020) are shown on the left, while the strains due to the second peak are shown on the right (September 2020). An equal strain response is observed near the niches due to excavation. The change of saturation induces larger deformation in the vicinity of the open niche. Due to anisotropic elastic properties, the strains in vertical direction are larger than in the horizontal direction

prototype developed by Graebling et al. (2022a) and briefly discussed in Sect. 2.6. As aforementioned, the bedding direction follows the horizontal direction of the niches' plane. With respect to the gallery, the bedding shows an average orientation of (147/57) degrees (Jaeggi et al. 2020; Ziefle et al. 2021) as pointed out in Fig. 2.21.

For a comprehensive and comparative visualization, we plot the saturation behavior in the model, Fig. 2.22. The upper part of the model is not shown in the visualization. Similar to the two-dimensional model, the saturation decreased in March 2020. A saturation decrease is observed in the open niche, lock and gallery as a result of the pressure changes. Due to the application of a constant pore pressure (boundary condition) in the closed niche, the saturation remains very high, i.e approximately 1. Note that the extent of the de-saturated zone in the vicinity of the gallery, lock and open niche remains approximately constant in size. An excavation damaged zone is not set explicitly in the model.

**Fig. 2.21** Setup of the three-dimensional model of the CD-A experiment. In-situ pore pressure and stresses are applied as initial conditions. The closed niche is subjected to a constant pore pressure, while the latter varies seasonally in the open niche

**Fig. 2.22** Saturation response in the three-dimensional model of the CD-A experiment. The saturation increases from March 2020 **a** to September 2020 **b** in the open niche and gallery

The saturation changes and a comparison between the modeled deformation and monitored convergence data is found in Ziefle et al. (2022c). As also observed in the 2D setting, the strain response of the model is affected by the excavation and by the seasonally change of air humidity. The strain response during the first de-saturation peak in March 2020 and re-saturation in September 2020 are shown in Fig. 2.23. Overall, all shown strain components increase up to winter and decrease in summer as observed, for example, in the strain component parallel to the axial direction of the gallery, Fig. 2.23a, b, and shear component, Fig. 2.23e, f. In the regions were the gallery and niches intersect, there is an increased influence of the saturation changes happening in the gallery. This can be visualized in all shown strain components, but specially in the direction parallel to the niches' axis as depicted in Fig. 2.23c, d. Note that there is a development of positive strains (stretching) on the walls of the niches and gallery. These contribute to an increase of the energy required for crack

**Fig. 2.23** Deformation response of the three-dimensional, anisotropic model. The results concerning the first desaturation peak (March 2020) are shown on the left, while the ones during the first de-saturation peak are shown on the right (September 2020). The strains are positive at the gallery and niches' walls. Note that the strains are larger in the open niche

propagation and contributes to the fact that cracks are prone to develop on the walls of the niches.

The described model behavior and mechanisms serve as a basis for model extension with respect to cracking. The proposed, extended model (Sect. 2.4.1) is first applied in a laboratory scale test related to the desiccation of a slab, which has been chosen due to the similarities with CD-A. Several tests concerning the energy splits, model formulation and understanding the influence of the numerical parameters involved in the coupled Richard mechanics phase-field approach have been carried out in Cajuhi et al. (2022). In the latter, the calibration and applicability of the coupled framework in the field scale are described in detail. The implemented framework is available in the open source software OpenGeoSys (Bilke et al. 2019, 2022). A general application of the framework in three-dimensional setting using a homogeneous geological window of the open niche and in context of desiccation cracking is discussed in Sect. 2.5.2.3. Before entering the fully coupled problem, we model and discuss the cracks initiated during the excavation using a FDEM approach in the following.

#### **2.5.2.2 Excavation Damaged Zone**

The Excavation Damaged Zone (EDZ) is considered as one of the main factors compromising the integrity in the near field of underground repositories. In the framework of FDEM, simulation of EDZ formation is captured with two separate model runs (Lisjak et al. 2014). In the first run, the in-situ stress state is applied by translating the specified stresses into nodal forces. After reaching the steady state, the displacement of outer boundaries are fixed allowing the model to maintain the in-situ stress conditions (Mahabadi et al. 2012). The excavation step would then be initiated employing a core replacement technique (Kavvadas 2005). Within this approach, the Young's modulus of the core material, i.e., the excavation region, is reduced in a step-wise manner to mimic the gradual de-confinement of the rock mass resulting from the tunnel advancement.

Due to the heterogeneous geology of the Opalinus clay rock, we evaluate the influence of the bedding layer orientation in forming an EDZ. In the first case, as shown in Fig. 2.24, the numerical model setup includes a 15 m × 15 m square domain with the opening region at its center with internal radius of 1.15 m. The geological window is representative of a one-quarter model of CD-A twins in which the bedding direction is horizontal. Additionally, the EDZ formation is studied with bedding direction inclined to tunnel axis at 45 degrees. In the second case, due to the asymmetric nature of the problem, a full 30 m × 30 m square domain has been considered. The employed material properties for transversely isotropic Opalinus clay is assumed same as the ones reported in Table 2.3. For both numerical test cases the vertical and horizontal in-situ stress are equal to 6.5MPa and 4.5MPa, respectively.

Figures 2.25 and 2.26 show the gradual increase of the failure pattern at different ratios of deconfinement, i.e., λ. Deconfinement is defined as the ratio between the decreased stress magnitude at the vicinity of the tunnel face to the initial magnitude of in-situ stress. Note that the crack propagation is parallel the bedding orientation in both cases. At the end of the simulation, we notice that the crack network formed with the horizontal orientation is more complex and connected in comparison to the case where the bedding is inclined with 45◦. On the other hand, the fracture network extends deeper into the formation for the second case. The percentage of broken joints as a function of distance from tunnel face are plotted for both cases at λ = 0.6 in Fig. 2.27a. The steeper curve represents the 45◦. Moreover, in Fig. 2.27b we show how the distribution of broken joints can be used in the creation of a model with

**Fig. 2.24** Two-dimensional FDEM geometry of CD-A experiment for modeling the excavation damaged zone

**Fig. 2.25** FDEM results of EDZ failure pattern near the open niche with horizontal bedding layer at **a** λ = 0.25, **b** λ = 0.75

**Fig. 2.26** FDEM results of EDZ failure pattern near the open niche with 45*circ* bedding layer at **a** λ = 0.3, **b** λ = 0.6

permeability values that depend on the distance from the niche wall, i.e., the near the wall the more permeable the model, since more joints are broken.

In other words, the complex crack network induced by the excavation can represent a zone of more instability and larger permeability in comparison to the intact rock. To further understand the latter, hydraulic models need to be taken into account

**Fig. 2.27 a** Distribution of broken joints percentage around underground opening. The percentage of broken joints in the model with horizontal bedding direction is larger (y-axis) than that in the 45◦ model. In the former, the cracks are less deep (x-axis) than those generated in the inclined model. Crack development motivates the EDZ abstraction concept and altered permeability profile shown in (**b**)

and coupled with the mechanical formulation. Ongoing work is being performed to account for hydraulic effects in the FDEM approach.

#### **2.5.2.3 Desiccation Cracks**

In this subsection, we use the fundamental equations from the continuum based methods (Sect. 2.4) to simulate desiccation cracks in three-dimensional setting. For the latter, no changes in the phase-field model are required, i.e. the same evolution equation is utilized and there is no need to adapt the algorithm to initiate the cracks or due to branching. The verification of the framework is based on the setup discussed in Cajuhi et al. (2018).

Since the open niche shows more pronounced hydro-mechanical effects, we propose a quarter three-dimensional model of a so called homogeneous geological window as shown in Fig. 2.28. The thickness of the model is the same as the radius of the niche (1.15 m). The total height/width is 5 m. A drying pore pressure related to the RH drop during the first winter is applied at the niche wall. At the remaining boundaries, no flow and symmetry boundary conditions are applied. Since the (free) movement of the niche wall, imposed by mechanical boundary conditions, can affect crack initiation and propagation and these conditions are not well defined in the field, we propose a free and a restrained setting. In the latter, the movement of the nodes at the front and back surfaces of the niche is hindered in the direction normal to the plane.

In the following, the results concerning the degree of saturation and cracking are shown. For visualization purposes, a threshold filter is applied on the crack phase-

field variable so that the elements where damage values equal or greater than 99% are removed from the blended out from the mesh. When the wall is allowed to move freely, a single crack initiates and propagates along as shown in Fig. 2.29a. The crack initiates at a pore pressure corresponding to approximately 94% RH and 80–81% saturation. As expected, restraining the front and back surfaces, and consequently the movement of the wall, induces earlier crack propagation and a more complex crack pattern. In the latter, the crack initiates approximately at 97% RH and 90– 92% saturation. Both cases lie in the expected drying RH range (Fig. 2.1), also reported in Regard et al. (2021), where crack development at the niche walls has been monitored and quantified.

The current phase-field evolution equation is not able to account for the mechanical and hydraulic anisotropies of the material. For this reason, the developed cracks do not necessarily propagate parallel to the expected bedding orientation as observed experimentally. Ongoing work is being carried out to extend the numerical approach with to anisotropic setting. Basis for the extension is the framework proposed in Ziaei-Rad et al. (2023) and discussed in Sect. 2.4. This extension should affect the crack orientation and induce cracks parallel to the bedding.

Additionally, the numerical simulation of desiccation cracking has been addressed in the framework of FDEM with coupled RM approach. The problem domain includes a 15 m × 15 m two-dimensional quarter model with 1.15 m niche radius. The assumed mechanical properties are the same as the ones considered in previous FDEM models dedicated to numerical simulation of fracturing in OPA. Fluid properties, porosity, values of permeability matrix and van Genuchten model parameters are all taken from Table 2.4. As mentioned in Sect. 2.4.3.2, the volume conservation equation and the mechanical equilibrium both are solved with finite difference explicit time integration scheme. This, consequently, forces smaller time steps for stability and

**Fig. 2.29** Joint visualization of saturation and crack phase-field responses of quarter model of the open niche. The model desaturates near the wall of the open niche (red color corresponds to drier conditions, while blue color to wetter conditions). For a direct visualization, the elements with damage equal or greater than 99% are not plotted, i.e. the empty spaces at the wall correspond to cracks. A single desiccation crack initiates at the wall of the open niche in the free setting. In the restrained setting, a more complex crack pattern develops

**Fig. 2.30** FDEM results of desiccation fracturing around the open niche after **a** 12 h and, **b** 24 h simulation time

convergence of the numerical solution. In the FDEM model, we impose a drying flux of *q* = 6.0 × 10−<sup>9</sup> m/s that acts on the niche wall for 24 h. The failure pattern as well as the saturation distribution after 24 h is shown in Fig. 2.30.

As can be seen, unsaturated areas, as a result of desiccation boundary condition, are only limited to the surrounding of the niche because of low permeability of Opalinus clay. Cracking starts from the top of the niche after approximately three hours of simulation and at approximately 3.5MPa suction (RH between 97–98%).

#### **2.6 Data Analytics and Visualization**

A prototype of the CD-A experiment using the 3D model presented in Sect. 2.5.2.1 has been prepared and discussed in detail in Graebling et al. (2022a). The model is placed in the Mont Terri rock laboratory and its results can be accessed via the proposed interactive tool. Furthermore, experimental data and its borehole positions can be combined and visualized together with the numerical model. In Graebling et al. (2022a), convergence measurements obtained from the boreholes via extensometers have been compared with the strains computed with the numerical model. The direct visualization and comparison of the data, combined with the experience to be able to understand where the data come from, i.e. by visualizing the borehole positions and the position of the experiment itself, increases the possibilities for exchange between different fields of expertise and with that, a multidisciplinary approach. Further details on the approach and examples of the visualization are discussed in Chap. 5.

#### **2.7 Summary, Conclusions and Research Transfer**

#### *2.7.1 Summary*

This chapter presented a workflow for understanding the Opalinus Clay characteristics, especially with respect to its shrinkage behavior. We carried out several laboratory experiments concerning the mechanical characterization of the rock. The results have shown and confirmed the strong influence of the anisotropy on the fracture behavior. In the SCB samples where the initial crack is parallel to the bedding, the crack propagates more easily, at a lower peak load, when compared to the experiments performed in the samples with initial crack perpendicular to the bedding. The tests have also shown the combined interaction of bedding and heterogeneities in the OPA since the cracks do not propagate uniformly in all samples. We also computed these experiments using FDEM and phase-field models coupled with the mechanics. Experience has also been gathered with respect to sample preparation and cutting procedures. For increasing the understanding on the effect of humidity variations, tests under different humidity conditions have been carried out. The tests were performed in a customized desiccator cell which has been equipped with different sensors to measure the controlled RH and temperature. In summary, the experimental tests show that an increase of suction (decrease of RH) leads to an increase of the strength of the material—as long as no fractures are induced that lead to an overall weakening of the material. The obtained data can be used as input in the simulations of drying/wetting and will affect the fracture behavior of the material in coupled hydro-mechanical simulations.

The laboratory tests have a direct connection with the field, since the samples were obtained from drilled boreholes in the vicinity of the in-situ CD-A experiment. In this region, the main lithology is characterized by the sandy facies of the OPA. We presented and used material properties obtained from the literature and from the experiments performed in laboratory. These data have been used as input in the various presented numerical models. The monitored data at the CD-A twin niches have been used for comparison with water content and convergence data. These results have been described briefly and are object of further investigations (Ziefle et al. 2022b, c). The twin niches at the CD-A experiment provide a good platform for comparing the hydro-mechanical effects under seasonal (open niche) and controlled high RH conditions (closed niche). Furthermore, the niches are equipped with several monitoring sensors, e.g. extensometers, permeability, water content, crackmeter, among others. The open niche shows more pronounced shrinkage effects in comparison to the closed niche, for example with respect to desiccation crack development.

The development of excavated damaged zones and further hydro-mechanical effects near the excavations, for example the initiation and propagation of desiccation cracks during ventilation, are motivating facts for developing and extending numerical approaches with respect to cracking and evaluating their applicability in the in-situ scale. With this objective in mind, we proposed continuum and discontinuum based models. A comparative study between these approaches for the laboratory scale is shown in Sattari et al. (2022). We described a continuum approach for modeling the unsaturated HM-coupled problem that has been implemented in Bilke et al. (2019, 2022). This approach has been applied in two- and three-dimensional settings. Aim of the simulations was to understand the effects of de-saturation and re-saturation in isotropic and anisotropic conditions for the hydraulic (permeability) and the mechanic (elasticity properties). The isotropic properties have been computed as an average of the anisotropic ones. We have shown the results for the first de-saturation peak, i.e. largest suction, and for the first re-saturation peak, i.e. largest pore pressure. The de-saturated zone increases up to the smallest pressure and part of it remains up to the re-saturation peak. We have observed the influence of the saturation changes on the deformation behavior and have stated that these are more pronounced in the open niche, where the seasonal behavior plays an important role and desaturation is stronger. These conclusions are valid in both isotropic and anisotropic models. However, in the latter, the desaturated zone extends further due to larger permeability along the horizontal direction.

The numerical results concerning the hydro-mechanical effects in 2D setting have shown that the contour and vicinity of the open niche are drier in comparison to the domain as a whole. This represents a driving mechanism for desiccation cracks. To obtain further information on the drying behavior of the walls, the model has been extended to 3D. This model is composed of open niche, closed niche, lock to the closed niche and gallery. We have observed an increase of the strains on the walls of the niches and gallery, with more pronounced effects in the latter and in the open niche.

Note that the models discussed up to now have not considered the development and extent of an EDZ. For modeling the crack networks formed during the excavation, we have used an FDEM approach. The approach combines the finite element analysis with a discontinuous representation of the crack, i.e. through zero-thickness joint elements that are inserted in the mesh. We observed a more complex crack network when modeling an horizontal bedding layer on the cross section of the niches that might represent a zone of weaker stability and larger permeability.

The unsaturated HM approach has been extended to account for cracks due to drying. The model has been applied into an homogeneous geological window, consisting of a quarter of the open niche. The free and restrained niche wall movement have been simulated to account for the full range of boundary conditions in the field. In the former, a single crack has been formed, while the latter induced a more complex crack pattern. The crack propagation is not direction dependency The effect of the hydraulic and mechanical anisotropies have not been taken into account in this model. The obtained crack patterns do not follow the bedding direction since the model is isotropic. Nevertheless, the quarter model has been able to predict desiccation cracks in the expected and monitored RH range (94–97%).

Finally, an interactive visualization of data in the field proportions has been proposed. The 3D model presented in this chapter has been inserted in this virtual prototype of the Mont Terri Laboratory, where the CD-A experiment takes place. The numerical results of saturation and deformation can be rendered and directly compared with borehole data, its position and monitored data provided by the sensors.

#### *2.7.2 Conclusions and Research Transfer*

We have proposed a workflow consisting of laboratory, field scale experiment and visualization of data. The laboratory experiments have provided new data related to the mechanical and fracture behaviors of the sandy facies of the OPA. Furthermore, a relationship between peak load and suction has been established for this material. The data obtained experimentally can be used as input for coupled hydro-mechanical simulations. For that, e.g. the strength needs to be determined for each tested suction range. An important model improvement would be to use strength-suction dependent material properties for simulating the material in laboratory and field scales.

The FDEM method showed good agreement with the measured laboratory experiments and has been also tested in the field scale for dealing with the EDZ. The numerically obtained EDZ shows complex crack networks, which might have an important effect on the permeability changes near the excavations. For further understanding this behavior, the hydraulic effects need to be coupled with the fracture-mechanical approach. This is object of ongoing work. A comparison between the response of the coupled hydro-mechanical methods (continuum and discontinuum based) is underway (Sattari et al. 2022). With that, we will evaluate the computational costs related to the inclusion of the discrete elements (joints) and of the addition of the phase-field evolution equation in each formulation.

The extension of the continuum framework based on the unsaturated HM approach using the phase-field model for brittle fracturing has the advantage and ability of simulating desiccation cracks. For a more resolved crack, the meshes need to be refined and this fact is directly connected to increased computational effort. Methods for automatic mesh refinement need to be further developed and explored in the context of HM coupled simulations. In the context of the unsaturated HM approach, the use of the presented anisotropic phase-field formulation is being evaluated and implemented in the open source framework. Furthermore, ongoing research is planned to account for swelling effects and crack closing as well as cracks due to swelling.

The long-term planning and execution of the CD-A experiment is an important aspect for testing the developed/extended modeling approaches and comparing its results to field data. Furthermore, the interactive data visualization of in-situ data is an important tool for communication and data exploration. The prototype presented in Graebling et al. (2022a) can be extended to account for further boreholes and sensors.

An interesting aspect for further research is the application of the presented methodology on other types of clays and rocks. This methodology can aid on finding the links between the different spatial and temporal scales. The combination of laboratory and field experiments, fundamental and applied research as well as interactive visualization represents an added value for different applications such as those in the geo-energy and radioactive waste disposal.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 3 Pathways Through Pressure-Driven Percolation in Salt Rock**

**Markus Knauth**

### **3.1 Hydro-Mechanical Properties of Salt Rocks and the Concept of Pressure-Driven Percolation**

Under undisturbed conditions and in the absence of significant impurities or intercalations (e.g. anhydrite) rock salts are considered to be tight against fluid or gas pressures, which stems from their crystalline structure composed of strongly viscoplastic, impermeable salt grains grown together at their grain boundaries (Minkley et al. 2013, 2015). Thus, its small natural porosity of 0.1–1% consists mainly of isolated brine pockets as evidenced e.g. by microscopic investigations of thin salt slivers. There are however stress- and pressure-conditions under which this hydraulically inactive flow network of connected grain boundaries can be opened. Then—and only then—a significant increase in macroscopic permeability can be observed due the creation of a connected flow path by opening the grain boundaries for permeation (Fig. 3.1). The conditions for such an increase in permeability can be classified into two principle mechanisms:


The first item describes the creation of flow paths by damaging the salt due to loading above its threshold for damage initiation (the so-called dilatancy boundary). The second mechanism—which is the main focal point of this study—describes the ability of fluids and gases to "open" the initially impermeable grain boundaries when the applied pressure becomes higher than the normal stresses acting on them. The attacking pressure thereby allows the medium to push its way into the grain boundaries and creates a connected flow path structure. Since this formation of a flow network occurs only after exceeding a certain threshold related to stress and pressure, this

M. Knauth (B)

IfG, Institut für Gebirgsmechanik GmbH, Leipzig, Germany e-mail: markus.knauth@ifg-leipzig.de

<sup>©</sup> The Author(s) 2023

O. Kolditz et al. (eds.), *GeomInt—Discontinuities in Geosystems From Lab to Field Scale*, SpringerBriefs in Earth System Sciences, https://doi.org/10.1007/978-3-031-26493-1\_3

**Fig. 3.1** Percolation of a colored solution on grain boundaries in hydro-mechanically coupled triaxial compression tests on salt specimen

effect has been termed "pressure-driven-percolation" in the field of salt mechanics. In the assessment of barrier integrity in salt mining—i.e. the preservation of aforementioned sealing capacity of the salt units—these two mechanisms have lead to the development of two associated criteria for their evaluation (Kock et al. 2012):


The dilatancy criterion describes the amount of mechanical damage the salt material has experienced, while the MSP compares the acting minimum principle stress with a theoretically attacking fluid pressure, e.g. a column of groundwater above the salt top or the gas pressure within a salt cavern. While the dilatancy criterion is typically more relevant for the assessment of damage at the contour of underground drifts and caverns (the excavation disturbed zone—EDZ), the MSP generally has a more far-reaching impact since the large-scale stress redistributions around a salt mine may reach far into the salt, potentially threatening its capacity as a hydraulic protection layer. It is important to note, that the minimum principle stress criterion is not the result of a hydro-mechanical modeling, but solely a stress- and pressure-based assessment of a given stress-state from a geomechanical modeling. As such it must be ensured, that it is applied at the most unfavorable state of the modeled system.

Additionally, the MSP is conservative in the sense that it only considers the scalar quantity of the least compressive stress. In reality, numerous laboratory hydromechanically coupled percolation studies have shown that the pressure-driven percolation is a highly oriented process, where the fluid generally percolates in a plane perpendicular to the least compressive stress (Kamlot 2009) (Fig. 3.2). Therefore it is possible that the MSP becomes too conservative, when the minimum principle stress may be lower than the attacking fluid pressure, but the direction of the stress field is such that this does not lead to a loss of the hydraulic protection capacity of the salt layer.

**Fig. 3.2** Highly anisotropic percolation in salt rocks perpendicular to the least compressive stress under true triaxial compression (Kamlot 2009)

Based on the previous considerations, the idea suggests itself to model the pressure-driven percolation as a diffusion process using a dynamic, full permeability tensor based on the fluid pressures and the stress tensor. However, this seemingly simple idea of modeling highly anisotropic diffusion leads to a number of numerical challenges in this and also in other scientific fields. The subsequent section will therefore outline the principle idea and the arising problems more detailed, before proposing a mesh-free numerical approach taken from the field of magneto-hydrodynamics in astrophysics.

#### **3.2 Development of a Mesh-Free Percolation Model Inspired by Magneto-Hydrodynamics**

#### *3.2.1 Numerical Challenges in Modeling Highly Anisotropic Permeation*

Most geomechanical codes capable of modeling a Darcy-type porous diffusion use finite-volume techniques (Moukalled et al. 2016), where the fluid resides in control volumes and each of these volumes may exchange fluid with its neighbors. By solving a balance equation and, the flow rates between those elements are determined and the corresponding changes in saturation and/or pressure are calculated. In order to determine the pressure gradient between to connected volumes i and j, a so-called two-point flux approximation (TPFA) is made (Aavatsmark 2002):

$$\nabla p\_{ij} = \frac{p\_j - p\_i}{d\_{ij}} \mathbf{e}\_{ij}$$

This widely-used approximation is one of the main reason for excessive numerical (i.e. false) diffusion when a highly anisotropic permeability tensor is used. For anisotropic permeability tensors, this approximation only yields accurate results when the directions of the element connection vectors (**e***i j*) are aligned with the Eigenvectors of that particular permeability tensor, i.e. its principle directions. Assuming that the permeability tensor is not constant and may change its principle directions (as it is proposed here by coupling it to the fluid pressure and stress tensor) this condition will never be met in the simulations. The introduced numerical error is of order *n*(0), meaning it is inherently tied to the model structure and does not converge away with increasing mesh resolution. It leads to numerical instabilities and implausible results including:


Even though this numerical problem is known for decades and the literature shows numerous correction approaches in a wide spectrum of applications from thermal propagation to astrophysics, it is still the default approximation in many—if not all established geomechanical codes. In the context of geomechanical barrier integrity not correcting for this problem leads to uncontrolled numerical diffusion, since the anisotropic fluid distribution will not stay stable and inevitably fill the whole model domain, thus rendering the approach useless for the evaluation of a hydraulic barrier.

#### *3.2.2 Basic Approach and Formulation*

In order to discretize the problem, a choice has to be made on how to discretize the model domain volume into a set of cells or particles i with coordinates *x<sup>i</sup>* . Instead of using a typical mesh e.g. of interconnected hexahedra or Voronoi elements, the approach formulated by Hopkins et al. (Hopkins 2015, 2016) considers a mesh-free alternative. A differential volume *d*<sup>3</sup> *x* at coordinate *x*. This differential volume can then be partitioned fractionally among its nearest particles by using a weighting function *W* depending on the distance *x* − *x<sup>i</sup>* and a "kernel size" *h*(*x*), thereby associating a fraction *<sup>i</sup>*(*x*) of the volume *d*<sup>3</sup> *x* with particle *i*:

$$\Psi\_i(\mathbf{x}) = \frac{1}{\omega(\mathbf{x})} W(\mathbf{x} - \mathbf{x}\_i, h(\mathbf{x}))$$

$$\omega(\mathbf{x}) = \sum\_j (\mathbf{x} - \mathbf{x}\_i, h(\mathbf{x})) $$

*W*(*x* − *xi*, *h*(*x*)) can be any continuous, symmetric and compactly supported function. Its absolute value is irrelevant due to the normalization. This type of model discretization therefore conceptually sits somewhere between a Voronoi tessellation and smooth-particle hydrodynamics. Hopkins et al. (Hopkins 2015) show, that this partition of unity discretization makes it possible to formulate the balance equation of property *U* for a given particle *i* with associated volume *Vi* as a Godunov-type finite-volume equation:

$$\frac{d}{dt}(V\_i U\_i) + \sum\_j F\_{ij} \cdot A\_{ij} = 0$$

with the "fluxes" *Fi j* in or out of an effective face area *Ai j* .

Another crucial advantage of this method is the fact that that type of discretization easily allows for a second-order accurate approximation of the gradient of any property *f* as a least-squares fit over all points in the kernel region, which effectively makes it a type of multi-point flux approximation method (MFPA):

$$(\nabla f)\_i^a = \sum\_j (f\_j - f\_i)(W\_i^{-1})^{a\beta}(\mathbf{x}\_j - \mathbf{x}\_i)^{\beta} \boldsymbol{\omega}\_j(\mathbf{x}\_i)$$

$$\mathbf{W}\_i^{a\beta} = \sum\_j (\mathbf{x}\_j - \mathbf{x}\_i)^a (\mathbf{x}\_j - \mathbf{x}\_i)^{\beta} \boldsymbol{\omega}\_j(\mathbf{x}\_i)$$

Is it important to note, that this gradient estimator is second-order accurate for an arbitrary particle configuration and therefore particularly suited for the mesh-free method.While the method proposed in Hopkins (2015) is constructed for very general hydrodynamics applications and leads to solving the associated Riemann problem, we now differ from the original formulation by transferring the general notion to the field of diffusion in porous media. In a simple explicit time-stepping scheme, we calculate the fluxes *Fi j* between each particle and its associated "neighbors" by evaluating the approximated gradient at the midway point between each particle and summing these contributions.

$$F\_{ij} = \frac{1}{\mu} (\mathbf{K} \cdot \nabla p\_{ij}) \cdot A\_{ij}$$

In order to stabilize the anisotropic diffusion equation, this flux is compared to the direct two-point flux

$$F\_{ij}^{direct} = \frac{1}{\mu} \left[ \mathbf{K} \cdot \left( \frac{p\_j - p\_i}{d\_{ij}} \mathbf{e}\_{ij} \right) \mathbf{A}\_{ij} \right],$$

Using this direct flux, the final flux *F*˜ *i j* is given by

$$
\tilde{F}\_{ij} = \begin{cases} 0 & \text{if } \text{sgn}(F\_{ij}^{direct} \cdot F\_{ij}) \text{ and } |F\_{ij}^{direct}| > 2 \cdot |F\_{ij}|, \\ F\_{ij} & \text{otherwise} \end{cases}
$$

Additionally, a basic saturation concept is introduced, giving each particle a porosity value and ensuring, that no flow can occur out of an unsaturated zone by multiplying the calculated flow rate with an empirical saturation function as *si j* used in the geomechanical modeling codes of Itasca (Itasca 2022a, b):

$$s\_{ij} = s^2 \cdot (3 - 2s)$$

where *s* is the saturation of the zone where the fluid is coming from (upwind scheme). Using this approach the volumetric flow rate *Q*˙*<sup>i</sup>* for particle *i* is

$$\dot{Q}\_i = \sum\_j \tilde{F}\_{ij} \cdot s\_{ij}$$

If the saturation at particle *i* is less than 1, this flow rate is used to update the saturation by

$$s\_i(t + \Delta t) = s\_i(t) + \frac{\cdot \mathcal{Q}\_i}{n\_i \cdot V\_i} \cdot \Delta t$$

with porosity *ni* and associated volume *Vi* . If the particle is already fully saturated, the associated pressure *pi* is updated instead using the fluid bulk modulus *Bf luid* :

$$p\_i(t + \Delta t) = p\_i(t) + \frac{B\_{fluid} \cdot \mathcal{Q}\_i}{n\_i \cdot V\_i} \cdot \Delta t$$

This algorithm was implemented as a Python module in FLAC3D, since this is the framework exposing the relevant internal variables of the corresponding numerical FLAC3D model in order to have access to the corresponding stresses and strains from the mechanical calculations of the salt behavior. In a staggered approach, the time-dependent mechanical problem is solved iteratively by FLAC3D for a given amount of time before the new stress tensor field is then given to the newly developed Python module for the anisotropic fluid flow calculation for the same timespan (albeit obviously using a different timestep for the fluid flow calculations). Although it is has experienced some optimizations, the Python implementation is—by nature of the interpreted programming language—much slower than a comparable C++ implementation using HPC-ready sparse matrix library would be. This would however need to be encapsulated again as a Python module in order to continue to work in conjunction with the FLAC3D software. Such efforts were tested in principle with the framework of this project and are (in addition to optimizations of the numerical structure itself) a potential field of improvement for future work. Together with the dynamic calculation of the local permeability tensor based on the stress tensor and the fluid pressure (see Sect. 3.3.3), we will subsequently refer to the new **M**esh-**F**ree **P**ercolation technique as the **MFP**-approach.

#### **3.3 Validation in Analytical Test Cases and Idealized Percolation Situations**

In order to demonstrate the capabilities of the new method and also show the unphysical numerical problems arising in standard implementations, a number of test cases have been assembled as a benchmark for the modeling of highly anisotropic porous flow. These examples and their comparative results are discussed hereafter and may aid in validation and verification for similar methods in the context of geomechanical percolation processes.

#### *3.3.1 Diffusion of a 1D-Step Function in Different Constant Permeability Fields*

In order to validate the basic diffusion capabilities, a 1D-step function is initialized in a constant isotropic permeability field (3.3) with a pressure *p* = 1 MPa on the left half and zero pressure on the right half of the model. The mesh has been created deliberately in a tilted fashion (with respect to the orientation of the quadratic elements) in order to create the most unfavorable condition for the subsequent analysis with an anisotropic permeability. For the isotropic permeability *k* in the first modeling, the pressure step will diffuse and soften the transition with time given by the analytical solution for the pressure profile *px*,*<sup>t</sup>* for a given time t:

$$p(\mathbf{x}, t > 0) = \frac{p}{2} + \frac{p}{2} \text{Erf}\left[\frac{\mathbf{x} - \mathbf{x}\_0}{4\kappa t}\right]$$

with

$$\kappa = \frac{\mathcal{B}\_{fluid} \cdot K}{n\mu}$$

where *n* is the rock porosity and μ the viscosity of the flowing medium. Figure 3.3 shows the pressure profile for a time of t = 600 s using the new MFP-approach in comparison with the native FLAC3D fluid logic. As to be expected, both codes are in good agreement with the analytical solution for the isotropic case.

In the second step, we now assign a anisotropic permeability, which is only nonzero in the vertical direction (*kx* = 0*m*2). Therefore, the numerical simulation should keep the initial step function perfectly intact, since there should be no flow into the x-direction. Again we examine the pressure profiles for both codes in Fig. 3.4 for a time of t = 12000 s. It can be clearly seen, that the MFP-approach coincides with the analytical (=initial) pressure profile, while the FLAC3D model exhibits the numerical errors discussed in the previous sections: The fluid front does not stay stable and move into the x-direction while also showing higher pressures than the initial condition on the left side and negative pressures on the right hand side. This is a clear example of the undesired effects using lower-order flow approximations

**Fig. 3.3** Model setup and mesh for the 1D isotropic diffusing step validation example (top) and the modeling results of the resulting pressure profile for different approaches (bottom) showing good agreement with the analytical solution

in highly anisotropic settings and is therefore proposed as a first and simple validity check for anisotropic approaches on porous flow.

#### *3.3.2 The Diffusing Ring Problem*

Another numerically challenging example for highly anisotropic flow is the so-called diffusing ring. Following the description in Hopkins et al., we initialize a purely

**Fig. 3.4** Resulting pressure profiles for the 1D anisotropic diffusing step function, which should stay perfectly immobile. The MHD approach successfully captures this, while the standard twopoint flux-approximation implemented in FLAC3D fails to keep the fluid front stable and also exhibits numerical problems

azimuthal pressure field in a 2D cubic box. With cylindrical coordinates centered on the origin we initialize *p*(*t* = 0) = *p*0(exp −(1/2)[(*r* − *r*0)<sup>2</sup>/δ*r* <sup>2</sup> <sup>0</sup> + φ<sup>2</sup>/δφ<sup>2</sup> <sup>0</sup> ]) with *p*0 =,*r*<sup>0</sup> = and φ<sup>0</sup> =. Since the permeability field is purely azimuthal, the pressure should then only diffuse along a ring-shape around the center. For short times, this problem has an exact solution:

$$p(t>0) = p\_0(\exp - (1/2)[(r - r\_0)^2/\delta r\_0^2 + \phi^2/\delta \phi^2])$$

with

$$
\delta \phi^2(r, t) = \delta \phi\_0^2 + 2Kr^{-2}t
$$

At later times, the diffusion from both sides of the ring intersects and there is no longer an exact solution. Again we test this problem setup using the MFP-approach and the native FLAC3D implementation. As an example of the modeling results, Fig. 3.5 shows the pressure distribution along the annulus after a short time (t = 60 s) in comparison to the original distribution. Even for this very short time, it can be seen that the MFP-results are again in good agreement with the analytical solution, while the FLAC3D results is already starting to deviate noticeably. For longer times, the pressure distribution wraps around as to be expected and stays fairly stable in the MFP modeling, while it again exhibits strong numerical diffusion in the FLAC3D implementation, eventually filling and equalizing the pressure in the whole domain. Therefore, this second—more elaborate—numerical example has confirmed

**Fig. 3.5** Resulting pressure profiles for the 1D anisotropic diffusing step function, which should stay perfectly immobile. The MHD approach successfully captures this, while the standard twopoint flux-approximation implemented in FLAC3D fails to keep the fluid front stable and also exhibits numerical problems

the shortcomings of the default methods and the capabilities of the MFP-approach to accurately model this challenging setup.

#### *3.3.3 Inclusion of Pressure- and Stress-Dependent Anisotropic Diffusion*

So far, the presented results were carried out for constant permeability problems with different boundary conditions and setups with the intent to verify the basic capabilities of the proposed method to accurately model anisotropic fluid flow. At this point we now include the mechanisms described in the introductory section in order to model the pressure-driven percolation in salt rocks.

In its principle axes, we formulate the anisotropic diffusivity tensor K in the following way:

$$K = \begin{bmatrix} k\_2 + k\_3 & 0 & 0 \\ 0 & k\_1 + k\_3 & 0 \\ 0 & 0 & k\_1 + k\_2 \end{bmatrix}$$

with

$$k\_i = \begin{cases} A(p + \sigma\_i) \text{ for } p > \sigma\_i\\ 0 \text{ otherwise} \end{cases}$$

where σ*<sup>i</sup>* are the principle stresses and *p* is the fluid pressure. By this construction, e.g. a pressure *p* exceeding the minimum principle stress σ<sup>3</sup> will lead to an increase

**Fig. 3.6** 45◦ degrees inclined pressure-driven percolation modeling with strong numerical diffusion and invalid pressures (holes) in the default FLAC3D implementation and stable percolation in the novel MHD approach

in permeability in the directions of σ<sup>1</sup> and σ2, i.e. the fluid may move in a plane perpendicular to the principle stress σ3.

Using this straightforward approach we can then test its functionality by modeling fluid injection tests in a cubic sample subjected to different stress conditions similar to the ones described in Kamlot (2009). On the purely cartesian mesh (20 × 20 × 20 elements), we initialize a fluid pressure of *p* in the center of the model, whose stress state is given by σ<sup>1</sup> = σ<sup>2</sup> = 11 MPa and σ<sup>3</sup> = 9 MPa. The given stress field was again deliberately chosen to lead to an unfavorably inclined fluid flow forming an angle of 45◦ against the principle axes of the mesh. Since the initial fluid pressure is below the minimum principle stress in the sample, there should be no flow occurring until the fluid pressure—which is increased stepwise—reaches and exceeds that mininum principle stress.

As expected, the initialized fluid stays stable until it exceeds the minimum principle stress. Upon exceeding the percolation front develops in the expected inclined orientation and is—most importantly—stable with respect to potential erroneous numerical diffusion in the direction of σ<sup>3</sup> (Fig. 3.6). Therefore, the proposed method is capable of both modeling both constant anisotropies as well as with the additional interaction of a spatially varying stress- and pressure-dependent permeability tensor.

#### *3.3.4 Application on a Large-Scale in Situ Borehole Pressurization Test in a Salt Mine*

A numerical recalculation of a large-scale in situ test was to be performed in order to demonstrate the capabilities and accuracy of the new approach. In this in-situ test, a 60 m long borehole with a diameter of 1,35 m was drilled between to mining

**Fig. 3.7** Location of the test site in the mining field Springen (left) and construction of the largescale 60 m borehole with 1.35 m diameter between mining levels (right)

levels of a salt mine and additionally monitored by a seismic array (Fig. 3.7). The borehole was cemented in the lower part, leaving a 40 m high cylindrical cavity which was then pressurized with gas. In accordance with the expected stress field, the salt borehole was gas-tight until the gas pressures reached a level above the minimum principle stress at the lower end of the borehole. The thereby created percolation front started in the bottom part of the borehole (at the cement plug) and progressed in a nearly horizontal plane above the underground mine drifts (Fig. 3.8), which was to be replicated by the newly developed method. The hydro-mechanically coupled calculations of the percolation reaction in response to the stepwise increase of borehole pressure successfully confirm the in situ observations. Using the novel MHD-approach, the modeled borehole is tight against the attacking gas pressure until the minimum principle stress is exceeded, which happens first in the lowest part of the borehole (Fig. 3.9). The direction of the percolation plane is oriented perpendicular to the smallest compressive stress and therefore initially nearly horizontally aligned. The percolation front in the in situ—test eventually hit a almost perfectly horizontal weakness plane and is therefore slightly more localized, which is not included in this modeling. However, the new method has shown to be well suited to model the pressure- and stress-tensor-dependent salt percolation process both in laboratory and in situ scale remarkably well.

#### 3 Pathways Through Pressure-Driven Percolation in Salt Rock 59

**Fig. 3.8** Pressure steps of the in situ—test (top left) and AE-measurement (setup top right) of the percolation front after the 68 bar pressure step in side (bottom left) and top view (bottom right)

**Fig. 3.9** Modeled percolation front starting to develop only after exceeding the necessary percolation pressure and permeating perpendicular to the least compressive stress

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 4 Pathways and Interfaces Under Stress Redistribution**

**Markus Barsch, Thomas Nagel, Holger Steeb, Patrick Schmidt, Dongwon Lee, Carlos Guevara Morel, and Jobst Maßmann**

There are several finite element-based modeling approaches to deal with fissures, fractures and discontinuities in rocks. A recent overview can found e.g in Mollaali et al. (2022). As one of them, lower-dimensional interface elements (LIE) were implemented in OpenGeoSys1 previously (Watanabe et al. 2012; Yoshioka et al. 2019) to map fractures and fissures discretely.

For example, in a 3D model consisting of a matrix and a fracture, the matrix is represented by 3D elements and the fracture by 2D elements. Crack growth can occur along pre-defined interfaces, i.e. the 2D elements, and therefore, the direction of crack path is predefined in contrast to other approaches resting on phase fields, extended finite elements or remeshing strategies.

1https://www.opengeosys.org.

M. Barsch (B) · T. Nagel

T. Nagel e-mail: thomas.nagel@ifgt.tu-freiberg.de

H. Steeb · P. Schmidt · D. Lee UoS, University of Stuttgart, Stuttgart, Germany e-mail: holger.steeb@mechbau.uni-stuttgart.de

P. Schmidt e-mail: patrick.schmidt@mechbau.uni-stuttgart.de

D. Lee e-mail: dongwon.lee@mechbau.uni-stuttgart.de

C. Guevara Morel · J. Maßmann BGR, Federal Institute for Geosciences and Natural Resources, Hanover, Germany e-mail: Carlos.GuevaraMorel@bgr.de

J. Maßmann e-mail: Jobst.Massmann@bgr.de

© The Author(s) 2023 O. Kolditz et al. (eds.), *GeomInt—Discontinuities in Geosystems From Lab to Field Scale*, SpringerBriefs in Earth System Sciences, https://doi.org/10.1007/978-3-031-26493-1\_4

TUBAF, Technische Universität Bergakademie Freiberg, Freiberg, Germany e-mail: Markus.Barsch1@ifgt.tu-freiberg.de

Further, lower-dimensional approaches taking into account the hydromechanics of deformable fluid-saturated fractures in porous media have been implemented in opensource software packages like DUNE2 or FEniCS3 (Schmidt and Steeb 2019; Schmidt et al. 2021). Numerically efficient coupling approaches for fractures and the surrounding porous rock allowing for HPC computations based on lower-dimensional models have been successfully applied to complex 2D/3D problems (Schmidt et al. 2022b) using modern code-coupling tools like preCICE.<sup>4</sup>

This general framework can be applied to a range of problems aside from fracture mechanics, including fracture slip problems and geotechnical problems with frictional interfaces, such as soil-structure interaction. In the sequel, we briefly summarize the theory and give a short insight into selected application examples from the GeomInt2 project. In cooperation with researchers from STIMTEC (Renner 2020) and others, those simulation tools have been used for the numerical interpretation of pressure transients obtained from large-scale experiments in the field, like e.g. harmonic pumping tests at the underground laboratory Reiche Zeche, Freiberg, Germany (Schmidt et al. 2021) or at the Grimsel Test Site, Switzerland, (Schmidt et al. 2022a). Based on those numerical-experimental investigations it was shown that fully-coupled hydromechanical approaches give a more physically sound insight into pressure transients especially if results are compared to simpler and more "classical" diffusion-based modeling approaches.

#### **4.1 Hydro-Mechanics of Deformable High-Aspect Ratio Fractures**

Viscous fluid flow in fractured porous rocks is causing strongly coupled hydromechanical interaction phenomena by local aperture, i.e. volume, changes of the fractures and resulting variations of the fluid pressure state. These phenomena have been investigated by various researchers over the last decades. For an overview we refer to a recent PhD thesis written during the period of the GeomInt and GeomInt2 project (Schmidt 2022).

Modelling approaches of fractured reservoirs require an efficient description of the hydro-mechanical processes. This includes on the one hand the accurate description of the deformation of the fracture and on the other hand the dynamic exchange of momentum between the fractures and the surrounding rock matrix. From a continuum perspective, model assumptions have to be made for the geometry of the fracture, and in a lower-dimensional setup, for the velocity profile within the fracture domain. For low Reynolds-number flow (Re < 1), a hybrid-dimensional model approach is obtained by an integration of the assumed parabolic flow velocity profile in normal fracture direction **e**ˆ3. This "a priori" integration is reducing the fracture flow domain

<sup>2</sup> https://www.dune-project.org.

<sup>3</sup> https://fenicsproject.org/.

<sup>4</sup> https://precice.org.

**Fig. 4.1** Representation of a single fracture -*Fr* embedded in a poro-elastic material body *<sup>B</sup>Pe* with material surface -*Pe*. For the kinematical description of the fracture, a local basis system **e**ˆ*<sup>i</sup>* is introduced. The basis vector **e**ˆ<sup>3</sup> is aligned with the fracture surface normal. The mechanical deformation of the fracture, and therefore the local volume change, is depicted with the evolution of the aperture δ, cf. Schmidt et al. (2022b), Schmidt (2022)

by one dimension, cf. Fig. 4.1. The resulting model is derived from 3D continuum mechanics by a consistent evaluation of the conservation of mass and balance of momentum under fully saturated conditions (Vinci et al. 2014; Schmidt 2022).

Without going into details of the derivations of models, we recapitulate here only the main governing equations. The bulk properties of the porous rock are described with the quasi-static linear poro-elastic equations (Biot 1941; Steeb and Renner 2019)

$$-\operatorname{div}\boldsymbol{\sigma} = \rho \,\mathbf{b},$$

$$\frac{1}{M}\frac{\partial p}{\partial t} - \frac{k^f}{\chi\_0^{fR}}\operatorname{div}\operatorname{grad} p = -\alpha \,\operatorname{div}\frac{\partial \mathbf{u}\_s}{\partial t}.\tag{4.1}$$

In Eq. 4.1 we introduced the local storativity 1/*M* of the porous rock, Biot's coupling parameter α, the isotropic hydraulic conductivity *k <sup>f</sup>* , and the effective weight of the fluid γ *f R* <sup>0</sup> . The total stresses are given by *σ*, the pore pressure is denoted with *p*, body forces are introduced as ρ **b** and the displacements of the solid skeleton are given by **u***s*. The set of coupled equations are supplemented with boundary conditions (Steeb and Renner 2019).

From first principles, a pressure diffusion equation could be derived for the evolution of the pressure state within the fracture. For high aspect ratio fractures, cf. Vinci et al. (2014), the final partial differential equation could be written in a simple non-linear format taking into account the evolution of aperture δ as well as fluid leakoff through an exchange term with the matrix *q*ˆ (normal component of the seepage velocity)

$$\frac{\partial \hat{p}}{\partial t} - \frac{\delta^2}{12 \,\eta^{fR} \,\beta^f} \text{div } \text{grad } \hat{p} + \frac{1}{\delta \,\beta^f} \frac{\partial \delta}{\partial t} = \frac{\hat{q}}{\delta \,\beta^f}. \tag{4.2}$$

In Eq. 4.2, we have introduced fluid properties like the effective dynamic viscosity η *f R*, the inverse of the fluid's bulk modulus β *<sup>f</sup>* and the fluid pressure in the fracture *p*ˆ. Details concerning the derivation and the mentioned dimensional analysis of Eq. 4.2 could be obtained from the references mentioned e.g. in Schmidt (2022). The set of coupled partial differential equations in weak form for the poro-elastic domain and the embedded fractures could be implemented into a Finite Element framework, e.g. within software packages mentioned in the introduction.

#### *4.1.1 Stiffness of Fractures*

The change of deformation of a fracture does not only depend on the fluid pressure transients *p*ˆ(**x**ˆ, *t*). The local deformation of the fracture is also affected by local contact forces between asperities of the solid surfaces. In case of natural fractures which have always rough fracture surfaces, those contact forces could for instance occur in the case of "mechanically closed" but still "hydraulically open" fractures. In the present lower-dimensional modelling approach, the corresponding effective fracture stiffness caused by local contact forces could be taken into account and expressed with the contact related normal stresses (Bandis et al. 1983)

$$
\sigma\_N^{\rm con} = E^{\rm Fr} \frac{1}{\delta\_{\rm max}^c - \delta^c} \delta^c. \tag{4.3}
$$

Note, that in Eq. 4.3 we considered only normal stresses, i.e. stress components normal to the fracture plane. Shear stresses within the fracture surface are neglected a priori. In the most simple linear constitutive approach, the effective contact normal stresses are proportional to the change of aperture described with the fracture closure δ*<sup>c</sup>* and the maximum fracture closure δ*<sup>c</sup>* max, cf. Gens et al. (1990) and results in Fig. 4.2 . In this linear constitutive relation we have introduced, as an inherent material parameter of the contact law, an effective normal stiffness *E*Fr.

On the laboratory scale, the introduced effective stiffness parameter could be determined in well-defined experiments. Therefore, natural fractures, obtained from mode 1 hydraulic fracture experiments, are mechanically tested, e.g. under harmonic excitation, cf. Fig. 4.3.

#### *4.1.2 Strong Versus Weak Coupling*

Implementation of the coupled set of partial differential equations introduced in terms of the biphasic poro-elastic and hybrid-dimensional flow model (Eqs. 4.1 and 4.2) might be distinguished regarding the chosen coupling scheme. In terms of the numerical strategy to reach equilibrium, chosen numerical coupling schemes are distinguished, cf. Fig. 4.4. Here, we may distinguish between staggered and

**Fig. 4.2** Left: Relative normal stress changes as a function of the effective fracture aperture. Normal stress changes are expressed relative to the acting normal equilibrium stress state which is equivalent to the sum of induced contact forces and the constant equilibrium normal stress state. Effective aperture fluctuations are induced by perturbations of the equilibrium state by means of fluid pressure changes or vice versa. Right: Consideration of the fracture's microstructure in terms of a contact factor, respectively controlling the difference between initial hydraulic and mechanical opening

**Fig. 4.3** Uniaxial setup for the determination of fracture stiffness *E*Fr under harmonic excitation at different amplitudes in combination with fracture roughness characterization through highresolution XRCT

monolithic approaches (Schmidt and Steeb 2019). Staggered numerical schemes treat the fracture and the rock bulk domain individually choosing solvers for both parts independently. Thus, highly efficient solvers for poro-elastic equations and for pressure diffusion equations could be adopted. This is especially interesting if large problems need to be solved, e.g. on HPC platforms (Schmidt et al. 2022b). Communication between the fracture domain *<sup>B</sup>Fr* and the the rock domain *<sup>B</sup>Pe* is conducted via boundary conditions, e.g. through special software packages like preCICE. Constraints have to be introduced in order to guarantee numerical stability.

**Fig. 4.4** Comparison of the weak coupling approach using non-conformal meshes (**a**) and a strong coupling approach with implemented interface elements (interface elements, and auxiliary nodes are explicitly shown for presentation purposes only; in the final discretization, the fracture surfaces align with *d* = 0) (**b**). Figure from Schmidt and Steeb (2019)

In contrast, monolithic coupling approaches build a single algebraic set of equations which is solved for both domains simultaneously.

Monolithic coupling of the fracture and the poro-elastic domain requires a mathematical discussion of the balance equations introduced to govern hydro-mechanical flow processes within the fracture, cf. Schmidt (2022), Schmidt and Steeb (2019).

#### **4.2 Lower-Dimensional Representation of Frictional Interfaces**

For constructing the weak form of the quasi-static equilibrium conditions

$$
\rho \operatorname{div} \boldsymbol{\sigma} + \varrho \mathbf{g} = \mathbf{0} \tag{4.4}
$$

in the presence of displacement jumps we introduce an enriched test function **v** consisting of a continuous (standard) part **v**<sup>c</sup> and a Heaviside enrichment *H*(**x**)**a** in the form

$$\mathbf{v} = \mathbf{v}\_{\mathbf{c}} + H(\mathbf{x})\mathbf{a}\_{\Gamma} \tag{4.5}$$

with

$$H(\mathbf{x}) = \begin{cases} -0.5 & \forall \mathbf{x} \in \Omega^- \\ 0.5 & \forall \mathbf{x} \in \Omega^+ \end{cases} \tag{4.6}$$

Therefore, (**v**<sup>+</sup> − **v**−)**<sup>x</sup>**- = [[**v**]]- = **a** represents the jump in the test function across the embedded discontinuity surface -. Superscripts + and − denote the opposite sides of the discontinuity surface. Then, we find the gradient of the test function as

$$\mathbf{grad}\,\mathbf{v} = \mathbf{grad}\,\mathbf{v}\_{\mathbf{c}} + \delta\_{\Gamma}(\mathbf{x}) \|\mathbf{v}\|\_{\Gamma} \otimes \mathbf{n}^{-}\tag{4.7}$$

where the Dirac function is defined and normalized as

$$\delta\_{\Gamma}(\mathbf{x}) = \begin{cases} 0 & \forall \mathbf{x} \notin \Gamma \\ \infty & \forall \mathbf{x} \in \Gamma \end{cases} \quad \text{with} \quad \int\_{\Omega} \delta\_{\Gamma}(\mathbf{x}) \, \mathrm{d}\Omega = 1 \tag{4.8}$$

To construct the weak form of Eq. (4.4)

$$0 = \int\_{\Omega} \mathbf{v} \cdot \left[ \operatorname{div} \boldsymbol{\sigma} + \varrho \mathbf{g} \right] \, \mathrm{d}\Omega \tag{4.9}$$

we perform partial integration of the stress divergence term by employing

$$\operatorname{div}\boldsymbol{\sigma} \cdot \mathbf{v} = \operatorname{div}(\boldsymbol{\sigma}\mathbf{v}) - \boldsymbol{\sigma} : \operatorname{grad}\mathbf{v}\_{\mathbf{c}} - \boldsymbol{\sigma} : \delta\_{\Gamma}(\mathbf{x})\|\mathbf{v}\|\_{\Gamma} \otimes \mathbf{n}^{-} \tag{4.10}$$

which results in

$$\int\_{\Omega \backslash \Gamma} \left[ \boldsymbol{\sigma} : \operatorname{grad} \mathbf{v} - \varrho \mathbf{b} \cdot \mathbf{v} \right] \, \mathrm{d}\Omega + \int\_{\Gamma} \boldsymbol{\sigma} \, \mathbf{n}^- \cdot \|\mathbf{v}\|\_{\Gamma} \, \mathrm{d}\Gamma = \int\_{\partial \Omega} \overline{\mathbf{t}} \cdot \mathbf{v} \, \mathrm{d}\Gamma \tag{4.11}$$

With *σ***n**<sup>−</sup> = **t** <sup>−</sup> and **t** <sup>+</sup> = **t**- = −**t** <sup>−</sup> we obtain

$$\int\_{\mathfrak{Q}/\Gamma} \left[ \boldsymbol{\sigma} : \operatorname{grad} \mathbf{v} - \varrho \mathbf{b} \cdot \mathbf{v} \right] \, \mathrm{d}\Omega - \int\_{\Gamma} \mathbf{t}\_{\Gamma} \cdot \|\mathbf{v}\|\_{\Gamma} \, \mathrm{d}\Gamma = \int\_{\partial \Omega} \overline{\mathbf{t}} \cdot \mathbf{v} \, \mathrm{d}\Gamma \tag{4.12}$$

Note that *σ* and **t** are the total stresses in a HM formulation of a fluid-saturated porous medium and thus couple into the following flow equations by means of the effective stress principle:

$$
\boldsymbol{\sigma} = \boldsymbol{\sigma}^{\prime} - \alpha\_{\text{B}} p \mathbf{I} \tag{4.13}
$$

$$\mathbf{t}\_{\Gamma} = \mathbf{t}'\_{\Gamma} - \alpha\_{\mathrm{B}}^{\mathrm{f}} p \mathbf{n}\_{\Gamma} \tag{4.14}$$

Note also that the present work assumes continuity of the pressure between matrix and fracture (*p*<sup>m</sup> ≡ *p*<sup>f</sup> ≡ *p*).

In the above *σ* is the effective stress tensor in the matrix, α<sup>B</sup> and α<sup>f</sup> <sup>B</sup> are coefficients weighting the pore-pressure contribution in the matrix and fracture to be defined suitably, *p* is the fluid pressure, **I** is the identity tensor, = φFR + (1 − φ)SR is the bulk density of the porous medium comprising the densities of the liquid and the solid, LR and SR respectively, of porosity φ and **g** is the gravitational acceleration. In the fracture with local unit normal **n**-, **t** represents the effective stress vector.

The effective quantities have to be supplied by constitutive equations. For the matrix, they are here given in incremental form

$$\text{d}\mathfrak{o}^{\prime} = \mathfrak{C}^{\text{m}} : \text{d}\mathfrak{e}^{\text{e}},\tag{4.15}$$

where *<sup>C</sup>*<sup>m</sup> is the material stiffness tensor, <sup>=</sup> <sup>1</sup> 2 grad **u** + grad <sup>T</sup>**u** is the linear strain tensor, <sup>e</sup> its elastic part and **u** is the solid displacement vector.

For the fracture effective stress, a connection is made constitutively to the displacement jump. Similar to the test function we specify **u** = **u**<sup>c</sup> + *H*(**x**)**w**(**x**-) such that (**u**<sup>+</sup> − **u**−)**<sup>x</sup>**- = [[**u**]]- = **w**. Then, in analogy to the incremental relationship for the matrix effective stress tensor, the fracture effective traction vector follows incrementally from elastic increments in the fracture face relative displacement vector

$$\mathbf{d}\mathbf{t}'\_{\Gamma} = \mathbf{K}^{\mathrm{I}} \mathbf{d}\mathbf{w}^{\mathrm{e}} \tag{4.16}$$

where d**w**<sup>e</sup> is the elastic part of the fracture relative displacement increment and **K**<sup>f</sup> is the second-order stiffness tensor of the fracture comprising normal and shear components.

For the following examples, the interface failure behaviour is of particular interest. Here, a Coulomb-type frictional interface was chosen

$$\mathbf{d}\mathbf{t'} = \mathbf{K}\mathbf{d} \|\mathbf{u}\| = \left[ K\_{\mathbf{n}} (\mathbf{n}\_{\Gamma} \otimes \mathbf{n}\_{\Gamma}) + K\_{\mathbf{s}} (\mathbf{I} - \mathbf{n}\_{\Gamma} \otimes \mathbf{n}\_{\Gamma}) \right] \mathbf{d}\mathbf{w}^{\varepsilon} \tag{4.17}$$

with a split of the displacement jump analogous to the strain split used in elastoplasticity **w** = **w**<sup>e</sup> + **w**p.

Based on the work by*Coulomb* in 1773 (Davis and Selvadurai 2002), the frictional failure criterion can be written as

$$
\sigma\_{\rm l} = c - \sigma\_{\rm n}^{\prime} \tan \varphi \tag{4.18}
$$

and describes a linear relationship between the limiting shear stress τ<sup>f</sup> in a plane and the normal effective stress σ <sup>n</sup> (tension positive) acting on that plane based on the cohesion *c* and friction angle ϕ of the fracture. The plane is in our case given by the fracture plane and a yield function *f*<sup>y</sup> describing the failure surface can be defined as follows

$$f\_{\mathbf{y}} = \tau - c + \sigma\_{\mathbf{n}}' \tan \varphi. \tag{4.19}$$

where τ is the acting shear stress in the fracture plane. For all states of stresses, which are within this envelope, i.e. *f*<sup>y</sup> < 0, no failure is predicted. This region grows as the fracture normal stress increases. Once *f*<sup>y</sup> = 0 is reached the fracture deformes plastically according to the flow rule

$$\mathbf{d}\mathbf{w}^{\mathbf{p}} = \mathbf{d}\lambda \frac{\partial g\_{\mathbf{y}}}{\partial \mathbf{t}'} \tag{4.20}$$

where the plastic potential is formulated to allow for non-associated flow:

$$\mathbf{g}\_{\mathbf{y}} = \boldsymbol{\tau} + \sigma\_{\mathbf{n}}^{\prime} \tan \psi. \tag{4.21}$$

**Fig. 4.5** Flow through discrete fracture network, where fractures are represented as 2D elements in a 3D setting (Lee 2022). Barsch & Nagel. Technische Universtität Bergakademie Freiberg. CC BY SA

where ψ is the dilatancy angle. Furthermore, the Kuhn-Tucker complementary conditions hold

$$f\_{\mathbf{y}}\lambda = 0 \qquad \lambda \ge 0 \qquad f\_{\mathbf{y}} \le 0 \tag{4.22}$$

From the implicit integration of the material model the consistent stiffness matrix is passed back to the finite element assembly routines.

#### **4.3 Verification and Application Examples**

The lower-dimensional interface approach can be used to simulate flow in fracture networks. An example from OGS simulations is shown in Fig. 4.5.

In this chapter, however, examples for mechanical and hydromechanical effects will be described briefly.

#### *4.3.1 Coulomb Model Verification*

To test the functionality of the Coulomb implementation, a compressive stress is created by compressing an element assembly separated by an interface vertically by a tenth of the initial fracture aperture, i.e. by 1 µm, in a displacement-controlled manner. The top element is then sheared horizontally into the positive *x* direction by 10 µm, followed by shear unloading and reloading into the negative *x* direction by 10 µm, while the vertical displacement is held constant. The latter boundary condition together with a positive dilatancy angle has the effect of increasing the fracture normal stress (accumulation of plastic normal fracture opening). Note that during this test the bottom boundary is fixed by a zero-displacement boundary condition whereas the boundaries at the left and right are freely movable.

Shearing causes a linear deformation in the horizontal direction until the applied shear exceeds the yield strength of the material. After this point, the shear stress increases at a lower rate until *ux* |*<sup>y</sup>*=2m reaches its maximum at maximum displacement, which illustrates the apparent hardening during the elasto-plastic deformation process under increasing confinement associated with the dilatant material response. Note that the material itself is not hardening but ideally plastic. Afterwards, the unloading /reloading due to shear reduces the shear stress in the same linear fashion observed in the initial stage of shearing (elastic unloading). After exceeding the yield strength, the material follows the dilatancy-induced apparent hardening again until maximum load reversal. The transition from linear deformation to material hardening after exceeding the yield strength clearly illustrates the effect of the Coulomb model on the overall deformation process. Also, the correct loading-unloading behaviour could be shown.

#### *4.3.2 Soil-Structure Interfaces*

Soil structure interactions play an important role in geotechnics. They include earth pressure development behind rough retaining walls, foundation tractions, interactions between tunnel liners and rock masses, driven piles, and many others.

Here, a simple laboratory test is chosen for illustrative purposes. Oedometer tests are a fundamental means of material characterization in geotechnical engineering. If the oedometer ring is compliant enough, radial stresses can be inferred in addition to the axial load in what is often called a soft oedometer. This example illustrates how the LIE model can be used to allow for relative displacement between the soil sample and the deformable walls of the soft oedometer.

The soft oedometer setup consists of two domains. The larger one represents the soil sample and the smaller one the compliant metal ring of the oedometer. The top boundary is assigned a constant normal traction representing the applied load. The bottom of the oedometer is fixed in the vertical direction and free to move in horizontal direction. The top and the right side of the metal ring are defined as traction free.

Due the Coulomb criterion, the soil can slide along the metal ring of the oedometer with a constant friction angle. Here, however, a frictionless interface was modelled. As illustration of the displacement discontinuity is the objective, the geomaterial is represented crudely oversimplified by a linear elastic model. The material properties of the benchmark are listed in the following table (Table 4.1):

Under the action of the top load, the soil gets compacted vertically (Fig. 4.6). The metal ring on the other hand experiences nearly no vertical displacements. The displacement jump in the solution is clearly visible by the transition from the smooth


**Table 4.1** Material properties

**Fig. 4.6** Vertical displacements in the soft oedometer (OGS simulation). The use of frictional lowerdimensional elements allows for the establishment of a vertical displacement jump between the soil (smooth color gradient, left) and the oedometer ring (red zone with zero vertical displacements, right) in the solution on a single domain

color gradient to the red outer part of the domain. This jump is inherent in the finite elements pace due to the Heaviside enrichment described above.

For the chosen parameters and numerical settings, the radial displacements causing the expansion of the ring used for radial stress measurements reach about 5.5 µm at the circumference and cause radial stresses of around 0.229MPa. This value is slightly lower than the 0.25MPa that would be expected for a rigid ring. Likewise, the axial compaction of the sample itself is somewhat larger than in the case of a rigid ring. The exact values can be confirmed analytically by solving the linear elastic equilibrium problem for the sample subject to the top load boundary condition in conjunction with the pressure vessel equation, where both domains are connected via radial stress and displacement continuity.

#### *4.3.3 Fault Slip Experiment (Mont Terri)*

Simulations related to the FS experiment in Mont Terri (Guglielmi et al. 2017) are described elsewhere (Rutqvist et al. 2020). Here, we show more recent results

using the LIE implementation described above. This example also links methods and application settings of WPs 2 and 3.

To simulate the fault slip displacement induced by fluid injection-induced overpressures, we define a 3D model domain (Fig. 4.7), which broadly represents the minor fault geometry embedded in the rock matrix (Rutqvist et al. 2020). The model domain has edge lengths of 20 m and contains the fault dipping at 65◦. Symmetry in the *x* direction is assumed and used to reduce the computational effort, therefore the side-length in the *x* direction reduces to 10 m.

The host rock matrix is modelled using linear elasticity assuming an isotropic material behavior with a bulk modulus of 5.9 GPa and a shear modulus of 2.3 GPa, based on the average values derived from experiments on Opalinus Clay at Mont Terri. For all outer boundaries of the matrix we assume roller boundary conditions. Based on a simplified experimental data set, the in-situ stress field is characterized by σ*zz* = −7MPa ,σ*yy* = −6MPa σ*xx* = −3.3MPa, where negative values indicate compression. Hydraulically, the host rock is considered impermeable because of the very low permeability of Opalinus Clay and the relatively short time-frame of these experiments. Fracture flow with various permeability models (constant, cubic law, cubic law only after shear slip, etc.) is applied while the fracture is assigned Coulomb behaviour. In this hydraulic-mechanical approach, a sudden increase in hydraulic aperture by a pre-set factor is modelled once shear failure occurs (Rutqvist et al. 2020). After fault activation, further aperture changes due to hydraulic-mechanical coupling are calculated by the cubic law. To parameterize the fluid phase and the porous medium we use standard values, see Table 4.2.

The initial fluid pressure estimated from site measurements was set to 0.5MPa while the effective fracture traction is calculated from the stress field in the matrix


**Table 4.2** Material properties for the modelling scenarios of the fault-slip experiment, scenario 1

**Fig. 4.8** Injection volume and fracture aperture change (numerical and experimental) using fullycoupled LIE models driven by step-wise injection pressure increases

by appropriate tensor rotation operations. The pressure and stress fields are assumed to be uniform and we neglect gravity in this case.

An initial comparison between the results of the numerical calculation and the laboratory data still shows an offset in the results (Fig. 4.8). Currently, parameter studies on the initial stress field, the hydraulic and mechanical fault properties, the activation mechanism and the injection control are undertaken to improve the results. Other modelling approaches were published in an overview paper (Rutqvist et al. 2020). Another modelling approach implemented in OGS was published in Urpi et al. (2020). These references contain details on the hydraulic models and the interpretation of the experimental data. A publication based on the approach presented here is currently under preparation.

**Fig. 4.9** Area of fault affected by shear slip

The area affected by shear slip, in which the fault transmissibility is drastically changed due to fault activation, is highlighted in Fig. 4.9. The area is centered around the injection point and is roughly circular in shape. This is inline with other modelling studies (Rutqvist et al. 2020).

#### **4.4 Concluding Remarks**

Lower-dimensional continua are mechanically and numerically intricate ways of representing interfaces of various kinds, including contact zones, fractures, faults, etc. Conceptually, the co-dimensional entities are generated by integration in their normal direction while making specific assumptions on the distribution of certain field quantities in this direction. Numerically, they can be represented by enrichment of finite element spaces and interface elements, for example. This kind of representation comes with both advantages and disadvantages compared to other methods. While this chapter showed examples on hydraulic, mechanical and hydromechanical behaviour of elastic and frictional interfaces, further aspects related to fracture mechanics and hydraulic fracturing are discussed in Yoshioka et al. (2019), Mollaali et al. (2022).

#### **References**


Yoshioka, K., Parisio, F., Naumov, D., Lu, R., Kolditz, O., & Nagel, T. (2019). Comparative verification of discrete and smeared numerical approaches for the simulation of hydraulic fracturing. *GEM - International Journal on Geomathematics, 10*(1), 13.

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 5 Virtual Reality and Computational Efficiency**

**Karsten Rink, Nico Graebling, Lars Bilke, Jörg Buchwald, Thomas Fischer, Christoph Lehmann, Tobias Meisel, Dmitri Naumov, Wenqing Wang, Keita Yoshioka, and Olaf Kolditz**

In this chapter we briefly describe information methods and technologies supporting geotechnical systems analyses, i.e. using virtual reality methods for data and model integration (Sect. 5.1) and improving computational efficiency by using highperformance-computing techniques (Sect. 5.2).

The generalised approach of the developed visualisation methodology makes it possible to integrate a variety of data formats into 3D scenes and to create studies for a wide range of application fields<sup>1</sup> For example, a prototype for an experiment information system was developed for the rock laboratory in Mont Terri operated by the Federal Office of Topography surveys Switzerland (swisstopo) (Bossart et al. 2017; Jaeggi et al. 2017). This prototype combines a representation of the complex geometry of the Swiss Jura mountains with tunnel system and boreholes of the rock

1http://www.ufz.de/vislab

L. Bilke e-mail: lars.bilke@ufz.de

J. Buchwald e-mail: joerg.buchwald@ufz.de

T. Fischer e-mail: thomas.fischer@ufz.de

C. Lehmann e-mail: christoph.lehmann@ufz.de

T. Meisel e-mail: Tobias.Meisel@web.de

D. Naumov e-mail: dmitri.naumov@ufz.de

© The Author(s) 2023 O. Kolditz et al. (eds.), *GeomInt—Discontinuities in Geosystems From Lab to Field Scale*, SpringerBriefs in Earth System Sciences, https://doi.org/10.1007/978-3-031-26493-1\_5

77

K. Rink (B) · L. Bilke · J. Buchwald · T. Fischer · C. Lehmann · T. Meisel · D. Naumov · W. Wang

UFZ, Helmholtz Centre for Environmental Research, OGS Core Team, Leipzig, Germany e-mail: karsten.rink@ufz.de

laboratory with a series of results of coupled numerical simulations for the planning and validation of long-time experiments, some of which are designed to run over several decades.

Due to the resolution, runtime and complexity of the simulations, the previously described methods for data reduction are essential for the visualisation of the results (Fig. 5.3). Figure 5.1 depicts a map of the URL Mont Terri including local labels for

W. Wang e-mail: wenqing.wang@ufz.de

N. Graebling · K. Yoshioka · O. Kolditz UFZ, Helmholtz Centre for Environmental Research, Leipzig, Germany e-mail: nico.graebling@ufz.de

K. Yoshioka e-mail: keita.yoshioka@ufz.de

O. Kolditz e-mail: olaf.kolditz@ufz.de

N. Graebling Leipzig University, Leipzig, Germany

D. Naumov TUBAF, Technische Universität Bergakademie Freiberg, Freiberg, Germany

K. Yoshioka University of Leoben (UoL), Leoben, Austria

O. Kolditz Technische Universität Dresden, Dresden, Germany the experiments which have been integrated into the Mont Terri visualization study. The development was carried out in close cooperation with swisstopo as well as the expert scientists who are in charge of the experiments.

The visualisation of geotechnical processes requires three-dimensional representation options for coupled field problems as well as validation in the context of corresponding measurement results, therefore, in the following we discuss aspects of data and model integration into a VR framework.

#### **5.1 Visual Data and Model Integration**

A prototype of a Digital Twin for URL Mont Terri has been designed and implemented (Graebling et al. 2022). In Sect. 2.6 the importance for data visualization and placing models into the real geological context has been already elaborated. In order to set-up the Digital Twin a particular Task VR (virtual reality) has been established in the portfolio of the Mont Terri project which was supported by the consortium. The VR task is dedicated to both data and model integration in order to combine the two sources of information for process characterization in the subsurface. The VR task has been evolved as a starting point of a new kind of data management for the future–a visually mediated data base including predictive capacities by integrated process models.

From the technical side, several software tools are important for the development of the VR study. The DataExplorer as part of OpenGeoSys-Workflows is used for data integration (Rink et al. 2022). The Unity framework is being adopted for interactive visualization purposes and data base access interfaces have been implemented.

#### *5.1.1 Data Integration*

An important prerequisite for data integration was linking two data bases, (i) the borehole information system BIS (developed a hosted by Simultec) and (ii) the sensor data base Geoscope (developed and hosted by Sixense Soldata). The support from the two companies in helping for the data base links is greatly acknowledged.

An illustration of data integration for URL Mont Terri is depicted in the figure collection. Figure 5.2a shows the entire tunnel system of URL Mont Terri including the positions of all active experiments (green spheres). The acronyms stand for the respective experiments. The indicated positions of the experiments can be clicked in the Unity framework and further information for connected boreholes and experiments become available. By clicking a particular borehole related geometrical information is displayed in a borehole menu (Fig 5.2b). Measurements such as displacements due to rock convergence can be mapped on the related boreholes of an experiment, therefore, structural and process data are being combined (Fig. 5.2c).

**Fig. 5.2** Data integration for URL Mont Terri: Prototype CD-A experiment

Structural information such as fracture planes from surface scans can be added into the scene (Fig. 5.2d). Other documentation such as related reports, publication links or photos from cores can be displayed in the context of the selected experiment or borehole (Fig. 5.2e). Moreover, Fig. 5.2f shows sensor data for the temporal evolution of displacements measured from different sensor positions. Therefore, the experimental information systems allows to show and combine structural and process data distributed in time and over space.

#### *5.1.2 Model Integration*

Model integration of three examples of Mont Terri experiments are presented:


Simulation results for the three experiments using OpenGeoSys are integrated into the VR study showing the process evolution in the real context on the URL Mont Terri.2 In case of the CD-A experiment, the numerical simulation has been used for the experimental design (calculation of optimal distance between the two niches) (Fig. 5.3).

#### **5.2 Computational Efficiency**

In order to achieve sufficient computational efficiency for the simulation of strongly coupled and often non-linear multi-field problems such as for thermo-hydromechanical (THM) processes, massive parallelization technologies have to be applied.

OpenGeoSys (OGS) is an open source software developed at the UFZ for the simulation of coupled thermal-hydrological-mechanical-chemical (THMC) multi-field processes in porous media. Such computations place high demands on memory and computing power. The aim is to model real-world applications by increasing spatial

<sup>2</sup> The Virtual Reality study for the Underground Research Laboratory Mont Terri (VR Task) has been supported by the Mont Terri consortium which is greatly acknowledged.

**Fig. 5.3** Model integration for selected experiments for URL Mont Terri

resolutions and to make models more realistic, for example by integrating non-linear material laws into the modelling. With a high spatial resolution, important geological structures that are relevant for the physical processes in the domain can be incorporated into the simulation model. This allows, for example, detailed simulations in large-scale areas as they are necessary for predictions about the mechanical integrity of various host rock types. The implementation required an efficient use of the available HPC resources. Within the framework of the project, scalability and parallel efficiency of individual OGS components (e.g. I/O, various local assemblers, monolithic approach vs. staggered scheme, linear equation system solvers) have been analysed. For the HPC optimisation, tracing/profiling tools were used to identify time-critical code sections and to find improvement potentials. Possible re-design steps for OGS were then derived from the analysis results and implemented. The optimised OGS code was successfully demonstrated in various application areas, such as hydrology and environmental geotechnics.

One of these time-critical code section was identified as the simulation output file operations. Before we used parallel file I/O based on the Visualization Toolkit (VTK) but its file formats were not designed with massive parallel clusters and file systems in mind and were a serious bottleneck in some application scenarios. We therefore implemented simulation result output based on HDF5 which decreased the time for writing simulation output files by more than an order of magnitude. Post-processing can still be done with common visualization tools like ParaView.

The variational phase–field model in this study has been implemented in parallel in OpenGeoSys. Unlike numerical methods that treat fracture explicitly with element edge or node splitting, the variational phase–field model treats fractures through a phase–field variable, which is a primary variable to be solved. Therefore, its parallel

**Fig. 5.4** THM simulations for a complete deep geological repository, shown is the temperature distribution 100 years after emplacement (Wang et al. 2021) (left), THM simulations for the fullemplacement (FE) experiment in the Mont Terri underground laboratory, shown are temperature, saturation and vertical stress distributions (Raith et al., 2020) (right)

implementation does not require any special consideration. If one has an implementation of linear elasticity, all that is required is to add an extra equation system for phase–field evolution and iterate the linear elastic system and the phase–field evolution system. All the variational phase–field simulations in this study were run in parallel and the linear scaling has been confirmed up to 100's of CPUs.

The application examples depicted in Fig. 5.4 present complex thermo-hydromechanical (THM) simulations for geotechnical problems such as the design of deep geological repositories for the storage of energy waste (e.g. heat-generating radioactive substances) in the geological subsurface. The HPC methodology is particularly demanding here because, on the one hand, it deals with non-linear coupled multi-field problems that place high demands on the parallel solution methods. Secondly, they involve complex geometries (mapping of geological and geotechnical structures), which place high demands on the domain discretisations and corresponding domain decompositions.

#### **5.3 Software Engineering and Benchmark Workflows**

#### *5.3.1 Containers*

Linux containers are a technology for encapsulating the runtime environment of software and executing it on any Linux system. Common technologies here are Docker and Singularity. Singularity is ideal for HPC environments because it enables seamless integration into parallelised environments (e.g. with MPI) and also uses a better execution concept than Docker (Singularity containers can be executed on the system with limited user rights whereas Docker containers require administrator

**Fig. 5.5** Schematic workflow for container generation

rights). Singularity is thus widely available on HPC systems, for example on the systems of the TU Dresden ("Taurus"), the Forschungszentrum Jülich ("Juwels") and the Helmholtz Centre for Environmental Research ("Eve").

In order to provide various software configurations of OpenGeoSys (e.g. different OGS versions, integration of third-party libraries to extend the range of functions) through Linux containers, a container creation framework "OpenGeoSys Container Maker" was implemented, which provides a user-friendly parameterisation of the container to be created. The framework is available to users via a web interface on the OpenGeoSys GitLab server (used for software version management, project management and software quality assurance). The container is created according to user specifications and made available as a download (Fig. 5.5).

The containers include the OpenGeoSys software in binary form as well as other tools helpful for the preparation and execution of simulations, such as computational grid generators or Python scripting interfaces. Singularity containers are regular files and can be transferred onto the HPC system easily after creation via the web interface mentioned above. Parallelised simulations can now be carried out with the help of the container. We used these container-based OGS workflows successfully during the project.

#### *5.3.2 Benchmarks and Jupyter Notebooks*

A new development within the project is the provision of container applications for Jupyter Notebooks.<sup>3</sup> In addition to the OGS core and external modules / libraries (e.g. MFront, PHREEQC, PETSc), these containers also contain the Jupyter Notebook server application and a number of Python packages which can be added as needed (Fig. 5.6). After starting the container, the Jupyter Notebook can be accessed as a browser application and OGS can be executed using notebooks. Jupyter Notebooks also form a new basis for OGS benchmark presentation and integration. New test cases are formulated and explained in Python-based Jupyter notebooks which can intermix script logic with explanatory text and images. Moreover, the large variety of existing Python tools can be used for pre- and postprocessing of OGS simulation

<sup>3</sup> https://www.opengeosys.org/docs/userguide/basics/jupyter-notebooks/.

results. Figure 5.8 shows the OGS benchmark gallery page4 which is organized according to the THMC process coupling hierarchy (Fig. 5.7).

#### *5.3.3 Benchmark Workflows*

The OpenGeoSys benchmark gallery (Fig. 5.8) is organized according to the THM/RTP process hierarchy, thermo-hydro-mechanical and reactive transport processes. A specific process class is represented by a tile showing a simulation result of the related process class. After clicking a tile, the available benchmarks of a process class comes visible and can be selected. Typically, an OGS benchmark starts with a short description of the problem and showing most important results for the benchmark test. All benchmarks are linked with the OGS project file (prj-file), therefore, the benchmark settings are directly available through the gallery. All benchmarks are part of the OGS quality assurance workflow which is continuously running all tests (benchmarks and so-called unit-tests for basic functionalities) after any code changes, automatically. For new benchmarks Jupyter notebooks are available for user convenience and user-specific pre- and postprocessing operations. New process

<sup>4</sup> https://www.opengeosys.org/docs/benchmarks/.

**Fig. 5.8** OGS Benchmark Gallery organized by process classes. Benchmarks of a specific process class are behind the tiles

classes and/or those with new benchmarks are highlighted as featured processes on top of the benchmark gallery.

Big Data applications cover a wide range of requirements for the analysis of complex data. This usually includes not only the analysis itself, but also the necessary pre-processing steps for data integration and preparation as well as methods for user interaction and evaluation. Therefore, it is necessary to create and execute complex workflows even for seemingly simple analyses. Modern HPC architectures are an ideal basis for providing customised and therefore high-performance working environments for different application requirements.

Large amounts of data not only reach the limits of computability, but above all of comprehensibility. Simulation processes are usually divided into three parts: in pre-processing, the simulation model is defined and input data for the simulation is created. Typically, the simulation as the second part generates large amounts of result data, which must then be processed and analysed in post-processing. Postprocessing is usually based on a visualisation of the result data and takes place either on front-end computers of the simulation cluster or on the user's desktop computer. Only in the latter case a transfer of the simulation data is necessary. The final analysed and processed data as a result of post-processing are far less extensive than the simulation results. Performing post-processing directly during the simulation (in-situ visualisation) makes it possible to bypass the communication/data transfer bottleneck and transfer only the already analysed data to the user's desktop computer.

The development, numerical realisation, validation and use of integrative modelling tools for the simulation of problem-specific, coupled processes as well as the efficient data presentation in the Big Data context make essential methodological contributions to the systematic analysis and prediction of processes in various areas of environmental science, especially with temporally and spatially large model dimensions. In the overarching conceptual and methodological approach to model and software development, the project results achieved by the grant recipient UFZ play a leading scientific role in the efficient and sustainable planning and management of the environmental science systems under investigation.

Work on the further development of the OpenGeoSys software platform was necessary in order to be able to reliably and efficiently carry out simulations of complex coupled processes of environmental science applications (hydrology, geotechnics) on large scales using HPC systems. This provided evidence that OpenGeoSys is suitable for modelling environmental science processes at real sites. In this context, the research work on parallelising the software should be mentioned as necessary and appropriate. With regard to scientific 3D visualisation, the work carried out on the development of workflows, numerical algorithms and software modules for the synoptic in-situ representation of heterogeneous model data from conceptually different sources has laid the foundations for the integrated consideration of structural information, results from measurement campaigns and results of numerical simulation, which were not available at the beginning of the project work.<sup>5</sup>

<sup>5</sup> OGS workflow development and High-Performance-Computing (HPC) capabilities have been supported by the ScaDS.AI and Pilot-Lab ESM initiatives which is greatly acknowledged.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

### **Chapter 6 Synthesis and Outlook**

#### **Olaf Kolditz, Tuanny Cajuhi, Ralf-Michael Günther, Holger Steeb, Frank Wuttke, Keita Yoshioka, Norbert Grunwald, and Thomas Nagel**

The principal interest of the GeomInt project consists of the investigation of effects on barrier integrity of three host rock formations: clay, salt and crystalline. The project focuses on distinct physical processes that can influence barrier integrity in these rocks, particularly those related to swelling and shrinkage, pressure-driven

O. Kolditz (B)

UFZ, Helmholtz Centre for Environmental Research, Leipzig and Technische Universität Dresden, Dresden, Germany e-mail: olaf.kolditz@ufz.de

T. Cajuhi BGR, Federal Institute for Geosciences and Natural Resources, Hanover, Germany e-mail: Tuanny.Cajuhi@bgr.de

R.-M. Günther IfG Institut für Gebirgsmechanik GmbH, Leipzig, Germany e-mail: ralf.guenther@ifg-leipzig.de

H. Steeb UoS, University of Stuttgart, Stuttgart, Germany

F. Wuttke CAU, Christian-Albrechts-Universität zu Kiel, Kiel, Germany e-mail: frank.wuttke@ifg.uni-kiel.de

K. Yoshioka UFZ, Helmholtz Centre for Environmental Research, Leipzig, Germany e-mail: keita.yoshioka@ufz.de

UoL, University of Leoben, Leoben, Austria

N. Grunwald UFZ, Helmholtz Centre for Environmental Research, Leipzig, Germany e-mail: norbert.grunwald@ufz.de

T. Nagel TUBAF, Technische Universität Bergakademie Freiberg, Freiberg, Germany e-mail: thomas.nagel@ifgt.tu-freiberg.de

© The Author(s) 2023 O. Kolditz et al. (eds.), *GeomInt—Discontinuities in Geosystems From Lab to Field Scale*, SpringerBriefs in Earth System Sciences, https://doi.org/10.1007/978-3-031-26493-1\_6

percolation and stress redistribution. In the first part of the project, methodologies were developed in both experimental and numerical areas (Kolditz et al. 2021b). Combined model-experiment studies (MEX) were used as the first synthesis tool. The problem classes (processes and rock types) were systematized and model and code comparisons were made for test examples (benchmarks) and laboratory experiments. A further result of the first project phase was also to systematically evaluate (advantages and disadvantages) the suitability of certain experiments and numerical methods for the individual or combined process classes (see above), which can influence the barrier effect, and to make further studies available (Kolditz et al. 2021a). In the second phase of the project, the methodical tools were used in particular for the analysis of laboratory and field experiments. Model chains (Sect. 6.1) and formationspecific workflows (Sect. 6.3) were developed as further synthesis tools, which are briefly presented below. Moreover, the GeomInt project has benefited from the outreach to various national and international activities such as the Mont Terri project, BenVaSim, DECOVALEX, and EURAD (Sect. 6.2).

#### **6.1 Mechanical Integrity of Clay Rocks**

#### *6.1.1 Model Chains*

Barrier integrity can be affected by a number of processes. In the case of gas permeation, relevant processes are highlighted in Fig. 6.1. Gases can arise e.g. from corrosion, degradation of organic compounds, setting of cement or pyrolysis. Gases initially spread out in dissolved form diffusely and advectively. At higher pressures, a gaseous phase can form, which propagates further as a multiphase flow. If the gas pressure continues to rise, mechanical damage to the rock with dilatancy and eventually propagating fractures can occur. The EURAD project is currently investigating whether this damage in clay rocks is reversible when the fluid pressures decrease and the plastic rock material converges and the fractures close again (the so-called selfsealing effect). This process chain shown in Fig. 6.1 can be mapped using new calculation methods (Grunwald et al. 2022). Dis- and continuum mechanical approaches, which were developed in GeomInt, are combined. The variational phase field (VPF) method is used for simulating fracture propagation processes (Yoshioka et al. 2022).

#### *6.1.2 Benchmarking*

For a reliable description of complex thermo-hydro-mechanical (THM) and reactive transport (RTP) processes, extensive tests of both the physical/chemical plausibility and the correct implementation in numerical simulators are essential. Figure 6.2 shows an sketch of benchmarks for THM/RTP models.

#### 6 Synthesis and Outlook 93

**Fig. 6.1** Model chain for simulating consecutive processes of gas propagation in clay rocks leading to damage barrier integrity functions (Marschall et al. 2005), extended by Grunwald

An extensive suite of test examples was created for the development of the OGS-TH2M model (Grunwald et al. 2022). A hierarchical concept was developed and implemented that systematically checks all (reasonable) process couplings based on the individual processes (Fig. 6.3). The extensive collection of benchmarks (*>*100 tested test examples) from OGS is directly available via the portal and offers users an ideal introduction to THM/CB modelling. Some of the OGS benchmarks are already available as Jupyter notebooks (Sect. 5.3.2) and can thus be integrated into other Python applications with the corresponding advantages of user-specific data analysis.

The OGS-TH2M model, a new process class, developed by Grunwald et al. (2022). A comprehensive benchmarking has been conducted with the BenVaSim (Lux et al. 2021; Pitz et al. 2022) and DECOVALEX projects (Sect. 6.2.1).

**Fig. 6.3** Benchmarking hierarchy for TH2M processes (Grunwald et al. 2022)

### **6.2 International Collaboration**

Fundamental developments from the GeomInt projects are widely used and deployed in various international collaboration activities in particular DECOVALEX, EURAD and the Mont Terri projects.

### *6.2.1 DECOVALEX*

DECOVALEX (DEvelopment of COupled models and their VALidation against EXperiments) is an international research and model comparison project, initiated in 1992, for advancing the understanding and modeling of coupled thermohydro-mechanical-chemical (THMC) processes in geological and geotechnical systems. DECOVALEX is an international community effort. OpenGeoSys is actively involved into all tasks of the current D2023 phase.<sup>1</sup> Task G is particularly dealing with the integrity of barrier rocks. This task is aiming to better understand reactivation of pre-existing discontinuities for brittle host rocks. In particular, the potential for existing features to undergo shear displacements and related changes in permeability as the result of coupled thermal, mechanical, hydraulic and chemical effects can all have significant impacts on repository safety functions (e.g., creating permeable pathways or, for very large displacements, mechanical damage of waste packages).

Task G is organised in four steps with increasing complexity (Fig. 6.4). The Task route starts from fracture mechanics (Step 1: M process) then diverting into hydro-mechanical (HM) and thermo-mechanical (TM) processes and finally combined into thermo-hydro-mechanical (THM) processes. These steps are described in the following including the concept, experimental design and data which will form

<sup>1</sup> https://decovalex.org/D-2023/overview.html.

the foundation for benchmarking exercises and experimental analysis. The research results are presented in the interim DECOVALEX report and the essence in Mollaali et al. (2022) which contains a comprehensive model and code comparison study including dis- and continuum mechanical approaches from all Task G members (nine international teams). Task G has been working on interpreting results from GeomInt laboratory experiments in Freiberg (Frühwirt et al. 2021) and Edinburgh (GREAT cell facility) (McDermott et al. 2018; Fraser-Harris et al. 2020).

#### *6.2.2 EURAD*

EURAD is a joint European project heading to safe radioactive waste management (RWM) through the development of a robust and sustained science, technology and knowledge management programme. It supports timely implementation of RWM activities and serves to foster mutual understanding and trust between Joint Programme participants.<sup>2</sup> GeomInt closely cooperated within EURAD particularly by contributing modeling know-how to work packages DONUT and GAS which are aiming at improving the understanding barrier integrity of clay formations in Europe suited as host rocks for deep geological repositories. GeomInt, therefore, actively helps in building European competences in clay systems understanding. Through EURAD national and international activities such as GeomInt, iCROSS, BenVaSim (Lux et al. 2021), DECOVALEX are bridged (Fig. 6.5).

#### **6.3 Workflows—Mont Terri Project**

One of GeomInt's main synthesis services is the development and implementation of workflows for the analysis of geotechnical systems. This is to be shown as an example for the joint work in the Mont Terri rock laboratory. Only the main features

<sup>2</sup> https://www.ejp-eurad.eu/.

**Fig. 6.5** International collaboration scheme—process-based cooperation for building European competences

of the workflow are shown here and the full description in Chap. 2 is referred to for details. (link Project geography)

The experimental basis is based on BGR field experiments in the Mont Terri rock laboratory and subsequent laboratory experiments. For example, drillcore samples were taken from the sandy facies of the OPA and analyzed in the geotechnical laboratory of the University of Kiel. The anisotropic, fracture and shrinkage behaviors were examined (observation/monitoring) in detail. The CD-A field experiment at Mont Terri has played a central role as a bridge between laboratory and field scales (ref). Much of the data from the URL Mont Terri was visually integrated into an information system as described in (ref) (Data Integration) The laboratory experiments (e.g. on the fracture behavior of the clay samples) were simulated using different model approaches. The essential processes and parameters that are relevant for the barrier integrity could be identified. (Process Simulation/Data Analytics) On this basis, a first transfer of the models to the field scale could be carried out. Exemplary, we were able to model cracks due to excavation and due to desiccation using discrete and diffuse approaches. In addition, an experimental design was carried out in advance for the planning of the CD-A experiment: the optimal distance between the two niches of the CD-A experiment was determined by numerical modeling and then implemented accordingly. This is an example of the use of models for the planning of complex and often expensive field experiments (Data Analytics/Knowledge Transfer). Finally, much of the experimental and modeled data and results have been integrated into the Visual Information System for Mont Terri and are available for

#### 6 Synthesis and Outlook 97

**Fig. 6.6** Workflow concept applied to clay rock systems of the rock research laboratory in Mont Terri

various planning, discussion and demonstration purposes (knowledge transfer) (Fig. 6.6).

In addition to the CD-A experiment, currently more experiments from the Mont Terri rock laboratory are integrated into the workflow concept, e.g. the FE (full scale emplacement heater) and FS (fault slip) experiments targeting at the "proof-ofconcept" for comprehensive workflows for entire research labs. The laboratory tests carried out in this phase of the GeomInt project, e.g. fluid percolation, will serve as basis for the application of numerical methods to model in the latter in-situ tests. Future work will be dedicated to further professionalize the workflow concept and transfer the know-how and software framework to other research labs.

#### **References**


and verification procedure. *Geomechanics and Geophysics for Geo-Energy and Geo-Resources, 8*(107).


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.